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ABSTRACT 

We study the evolution of planetesimals in evolved gaseous disks, which orbit a solar-mass star and 
harbor a Jupiter-mass planet at ~ 5 AU. The gas dynamics is modeled with a three-dimensional 
hydrodynamics code that employes nested-grids and achieves a resolution of one Jupiter’s radius in the 
circumplanetary disk. The code models solids as individual particles. Planetesimals are subjected to 
gravitational forces by the star and the planet, drag force by the gas, disruption via ram pressure, and 
mass loss through ablation. The mass evolution of solids is calculated self-consistently with their tem¬ 
perature, velocity, and position. We consider icy and icy/rocky bodies of radius 0.1-100 km, initially 
deployed on orbits around the star within a few Hill radii (Ru) of the planet’s orbit. Planetesimals 
are scattered inward, outward, and toward disk regions of radius r ^ ap. Scattering can relocate 
significant amounts of solids, provided that regions |r —a^j ^ 3i?H are replenished with planetesimals. 
Scattered bodies can be temporarily captured on planetocentric orbits. Ablation consumes nearly all 
solids at gas temperatures > 220 K. Super-keplerian rotation around and beyond the outer edge of 
the gas gap can segregate < 0.1km bodies, producing solid gap edges at size-dependent radial loca¬ 
tions. Capture, break-up, and ablation of solids result in a dust-laden circumplanetary disk with low 
surface densities of km-size planetesimals, implying relatively long timescales for satellite formation. 
After a giant planet acquires most of its mass, accretion of solids is unlikely to alter significantly its 
heavy-element content. The luminosity generated by solids’ accretion can be of a similar order of 
magnitude to the contraction luminosity. 

Keywords: accretion, accretion disks — hydrodynamics — methods: numerical — planet-disk inter¬ 
actions —planets and satellites: formation — protoplanetary disks 


1. INTRODUCTION 

Planetesimal accretion is a key process in the forma¬ 
tion of a giant planet. In the Core-Nucleated Accre¬ 
tion (CNA) scenario (Pollack et al. 1996), planetesimal 
accretion accounts for the formation of the initial core 
(Phase 1). Additional planetesimal capture during the 
slow accretion of the gaseous envelope (Phase 2) releases 
gravitational energy, which must be radiated so that the 
envelope can contract. As the captured planetesimals 
pass through the envelope and ablate, they leave be¬ 
hind solid grains which affect the opacity. Both the 
energy release by sinking solids and the opacity affect 
the planet luminosity and help determine the duration 
of Phase 2 (Hubickyj et al. 2005; Movshovitz & Podolak 
2008; Movshovitz et al. 2010). This ablated material 
will also affect the subsequent composition in the enve¬ 
lope (laroslavitz & Podolak 2007; Mousis et al. 2014). 
After the rapid gas accretion and ensuing contraction 
(Phase 3), when the planet has acquired most of its mass, 
additional accretion of solids occurs both onto the planet 
and in the subdisk surrounding the planet (the circum¬ 
planetary disk). This additional accretion adds relatively 
little to the planet itself, but can have important conse¬ 
quences for the formation of the regular satellites (e.g., 
Canup & Ward 2009; Estrada et al. 2009) and, possibly, 
for the occurrence of the irregular satellites. 
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As computing power has increased, studies of planetes¬ 
imal accretion have become more detailed. Pollack et al. 
(1996) calculated the accretion rate assuming that the 
planetesimals were uniformly distributed over the feed¬ 
ing zone. Inaba et al. (2003) used a statistical model, 
which simulated such effects as scattering and collisions 
between planetesimals. Modifications and improvements 
of this idea have been employed over the past few years 
(Kobayashi et al. 2010; Bromley & Kenyon 2011). More 
recently, D’Angelo et al. (2014) revisited the problem 
combining detailed calculations for the evolution of a 
swarm of planetesimals and for the structure of the 
planet’s envelope, up to the beginning of Phase 2. 

The last stage of accretion, after the planet has un¬ 
dergone its rapid contraction, presents special difficul¬ 
ties. The planet is massive enough to open a gap in the 
gas density distribution of the circumstellar disk, and 
this constrains the motion of the gas flowing toward the 
planet (e.g., Lissauer et al. 2009, and references therein). 
The inflowing gas and the associated circumplanetary 
disk, in turn, affect how planetesimals are delivered to 
and captured by the subdisk and by the planet itself. Ad¬ 
ditionally, planet-induced perturbations on the gas may 
impact the redistribution of planetesimals in the circum¬ 
stellar disk. In this work, we present the results of cal¬ 
culations that combine the three-dimensional (3D) gas 
dynamics with an N-body module to study these pro¬ 
cesses in detail, including the determination of temper¬ 
ature, ablation, and fragmentation of planetesimals. We 
consider the case of a giant planet of one Jupiter’s mass 
and model solids in the (initial) size range 0.1-100 km. 
Results are presented also for smaller bodies, down to 
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Table 1 

List of Symbols 


Symbol 

Definition 

{r, e, (j)} 

Stellocentric spherical polar coordinates 

'^mn,mx 

Min/Max grid radius 

^mn,mx 

Min/Max grid co-latitude angle 


Stellar mass 


Frame rotation rate 


Planet rotation rate 

Mp 

Planet mass 

Rp 

Planet radius 

Ru 

Planet’s Hill radius 

CLp 

Planet semi-major axis 

rp 

Planet position 

V, 

Gas velocity 

pg 

Gas volume density 

^g 

Gas surface density 

Pg 

Gas pressure 

^g 

Gas sound speed 


Turbulence kinematic viscosity 

ag 

Turbulence viscosity parameter 

{po, So} 

Circumstellar disk reference densities 

H 

Circumstellar disk pressure scale-height 

Tn 

Circumstellar disk temperature 

Tg 

Gas temperature 

Pg 

Gas mean molecular weight 

Pg 

Gas molecular dynamical viscosity 

{Lr, Lq^ L(jf\ 

Particle specific linear/angular momenta 

{rs, Vs,aA 

Particle position, velocity, and acceleration 

{Fd,^d} 

Aerodynamic drag force and acceleration 

Cd 

Drag coefficient 

M 

Mach number 

7^ 

Reynolds number 

B 

Biot number 

Ps 

Particle density 

Ms 

Particle mass 

Rs 

Particle radius 

es 

Particle emissivity 

Ls 

Particle specific vaporization energy 

Cs 

Particle specific heat 

Ts 

Particle temperature 

Xs 

Particle thermal conductivity 

Ss 

Particle isothermal depth 

Ps 

Particle mean molecular weight 

(Js 

Nominal material compressive strength 

Pv 

Particle vapor pressure 

Ter 

Particle critical temperature 

Pdy 

Dynamical pressure 

Rdy 

Particle break-up radius 


Rock volume fraction of mixed medium 


Conductivity efficiency factors of mixed medium 

{x,y, z} 

Planetocentric cartesian coordinates 

r 

Distance from the planet 

Te 

Circumplanetary disk effective temperature 

H 

Circumplanetary disk local thickness 


Rosseland mean opacity 


1 cm in radius. 

In what follows, we describe how the thermodynamical 
evolution of the disk’s gas and of the solids is calculated 
in Sections 2 and 3, respectively. The numerical pro¬ 
cedures are outlined in Section 4. Results for the disk 
evolution are presented in Section 5, and those for the 
evolution of planetesimals are presented in Sections 6 
and 7. We conclude discussing our findings in Section 8. 
Finally, further details and numerical tests are given in 
Appendix A and B. 

2. THERMODYNAMICAL EVOLUTION OE THE 

DISK 

We work in a reference frame whose origin is fixed to 
the star and which rotates about the origin at a rate U/, 


equal to the angular velocity of the planet around the 
star, ftp. Eor a planet on a circular orbit, ftp is equal to 
the mean-motion: 


ftp — 


G(M^ + Mp 


a 


3 

p 


(1) 


where and Mp are the star’s and planet’s mass, re¬ 
spectively, and ap is the planet’s semi-major axis. The 
planet-to-star mass ratio is Mp/M^ = 9.8 x 10“^. Table 1 
contains a list of the main symbols used in this paper. 

The circumstellar disk is represented by a spherical 
sector with an inner hole. Consider a spherical polar 
coordinate system {O; r, 0}, where r indicates the po¬ 
lar distance from the origin, O, the angle 0 is the co¬ 
latitude (6> = 0 is the north pole, 0 = 7r/2 is the mid¬ 
plane, and 7r/2 — 6> is the latitude), and the angle 0 is 
the azimuth. The disk volume is within the range given 
by ["^mn^^mx] ^ [^mn? ^mx] ^ ^TT, where T^nn — 0.4: Cip aud 
^mx = 4ap. We assume that the planet’s orbit lies in 
the disk’s equatorial plane, 0 = 7r/2, and that the disk is 
symmetric with respect to this plane. Consequently, only 
half of the disk volume needs to be simulated. Therefore, 
we set Omn — 27r/5 and = '7r/2. 

The disk’s gas is approximated to a viscous fluid of 
constant kinematic viscosity z/^, volume density and 
velocity Vg. In the following, all gas-related quantities 
will bear the subscript 'g\ The viscosity I'g is typi¬ 
cally assumed to arise from turbulence (of unspecified 
origin) within the gas, and needs not to be confused 
with the molecular viscosity, introduced below, which is 
much smaller in magnitude. We set Ug = 10~^ ap ftp, 
which corresponds to a turbulence parameter (Shakura 
& Syunyaev 1973) ag = 0.004 for our choice of the disk 
thickness. 

We generally assume that the circumstellar gas is lo¬ 
cally isothermal (the temperature depends only on r) and 
that the pressure is 

Pg = c.gpg. ( 2 ) 


The gas sound speed is Cg = {H/r)v}^, and vy^ is the lo¬ 
cal Keplerian velocity. The relative thickness of the disk 
above the equatorial plane, H/r, is taken to be constant 
and equal to 0.05. Therefore, the simulated disk vol¬ 
ume extends in the vertical direction for more than 5.2 
scale-heights, H. By assuming the equation of state for 
an ideal gas, the temperature in the circumstellar disk 
becomes 


/ pgmu\ 

V fcB J 


"9’ 


(3) 


where pg is the mean molecular weight of the gas, mu is 
the hydrogen mass, and ku is the Boltzmann constant. 
Since oc 1/r, the gas temperature in the circumstellar 
disk is proportional to 1/r as well. 

Given the thermal state of the gas, the disk region in 
which the gravity of the planet dominates over that of 
the star has a linear size on the order of the Hill ra¬ 
dius, Ru = ap[Mp/{3M^)]^^^, which strictly speaking 
represents the distance of the Lagrange point Li from 
the planet (to leading order in Mp/M^, e.g., Kopal 1978; 
Murray & Dermott 2000). In absence of gas (i.e., neglect¬ 
ing pressure and viscosity effects), this region is a solid 
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Table 2 

Disk’s Gas Constants 


star and the planet, from non-inertial forces, and from 
the drag force exerted by the gas 


H/r 

Vg- 


Ig 



0.05 

10-^ 

2.39 

1.4 

10-100 

10-12_10-11 


^ In units of alClp. 

^ Unperturbed and pg, in cgs units, at 5.2 AU. 


of revolution (around the planet-star axis) whose volume 
is only about a third of that occupied by the Hill sphere. 
Thus, the effective (volumetric mean) radius of the re¬ 
gion is 2Rii/3 (Kopal 1959; Paczyhski 1971; Kopal 
1978; Eggleton 1983). 

Table 2 summarizes the disk’s gas parameters, assum¬ 
ing = Mq. The disk’s reference densities at 5.2 AU 
are derived from evolution models discussed in Section 5. 
The mid-plane temperature of the disk at 5.2 AU, for 
the choice of parameters in Table 2 , is ~ 120 K. As ex¬ 
plained in Section 5.1, the gas temperature distribution 
is given by Equation (3) but, in the restricted region of 
the circumplanetary disk, it is modified according to sim¬ 
ple arguments based on local viscous heating, black-body 
heating by background radiation, and radiative cooling. 

Orbital migration of the planet is neglected. At the 
higher Eq considered here, the planet would drift inward 
at a speed ^ (D’Angelo & Lubow 2008). At the 

lower Eo, migration would be inertia-limited and slower. 

3. THERMODYNAMICAL EVOLUTION OF 
PLANETESIMALS 

Here we describe the physical model for the thermody¬ 
namical evolution of the planetesimals. Since the model 
generally applies to any solid particle, regardless of the 
size, in this section we shall refer to the planetesimals 
simply as particles. All particle-related quantities will 
bear the subscript ‘s’. Sometimes, for ease of notation, 
this subscript is dropped. 


3.1. Particle Dynamics 

Let us introduce the linear momentum per unit mass 
in the radial direction Lp = the meridional angular 
momentum per unit mass Lq = wq^ and the azimuthal 
angular momentum per unit mass = rsinOv^^ all 

defined in an inertial frame of reference. The velocity 
is the absolute azimuthal velocity: = T 0 +(r sin 6 >)^Uj. 

In terms of these momenta, the equations of motion of a 
particle can be written as 


dLp 

dt 


— ap -\— 

r 


dLo cosO ( L(k 

dt sm(J \rsm(J 

dL(j) . 

—= afhV smu, 

dt ^ 


r sin^ 
2 


( 4 ) 

( 5 ) 

(6) 


where a^, a^, and a(f) are the spherical components of 
the gravitational acceleration imparted to the particle. 
Notice that the subscript ‘s’ associated to the coordinates 
and momenta of the particle is dropped. 

In our case, the acceleration in Equations (4), (5), and 
( 6 ) arises from the gravitational forces exerted by the 




GMp{vp-Vs) GM, , 

|rp-r,|3 + 

GMp GMs 

— PtfX {PlfXVs) — 2flfXVs 


( 7 ) 


where is the position vector of the planet, r^, v^, and 
Ms are position, velocity, and mass of the particle. The 
rotation rate vector, Qf, is parallel to the direction of the 
north pole {0 = 0), i.e., flf = ftfZ. The components a^, 
a^, and are found by projecting along the spherical 
polar unit vectors. 

In Equation (7), the third term on the right-hand side, 
ai), is the drag acceleration. The fourth and fifth terms 
are the non-inertial accelerations imparted to the star 
(the origin) by the planet and the particle, respectively. 
The last two terms are the centrifugal and Coriolis ac¬ 
celerations. Additional terms may be included, such as 
the gravitational force per unit mass exerted by the disk 
on the particle, which we ignore here, and that exerted 
on the star, another non-inertial term, which is ignored 
as well. 

Notice that Equations (4), (5), and ( 6 ) use absolute 
linear and angular momenta, hence they apply regardless 
of whether the vector ftf is constant or not. Equation (7) 

is valid for Clf = 0 and requires the additional non- 

inertial term —CtfXVs in case ilf 7 ^ 0 . 

Let us indicate with As the cross section of a parti¬ 
cle, then the drag force experienced by the particle while 
moving through the gas is 


F 



CoAspglVg 


Vs|(Vg 


(8) 


where Gd i^ the drag coefficient. For spherical particles 
of uniform density ps and radius the drag force per 
unit mass becomes 


■*” = 5^ (ft) 

In general, the drag coefficient, Gd^ depends on the 
(relative) Mach number 


M = ( 10 ) 

^9 

and on the (relative) Reynolds number, which can be 
written as 

= ( 11 ) 

% 

in which pg represents the molecular dynamical viscos¬ 
ity of the gas. We use the expression for Gd derived by 
Melosh & Goldin (2008), which is a continuos function 
applicable over the full range of M. and IZ. In the contin¬ 
uum flow limit, that occurs when MIlZ 1 , we embed 
in the coefficient of Melosh & Goldin (2008) a drag for¬ 
mula proposed by Brown & Lawler (2003). The drag 
coefficient is discussed in more detail in Appendix A. 
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3.2. Particle Thermodynamics 

A particle moving through gas sweeps a mass per unit 
time equal to Agpg\vg — v^j. Collisions between gas 
atoms (and/or molecules) and the particle transfer some 
amount of the specific kinetic energy of the gas, \vg — 
V 5 P/ 2 , to the particle at a rate (/ek/ 2 )A^p^|v^ — 
where /ek is the fraction of the total collisional kinetic 
energy transmitted as heat to the particle. This energy 
exchange can also be interpreted in terms of the rate 
at which work is done on the particle by drag in the 
gas frame Fe)-(v^ - v^) = {CD/‘^)Aspg\vg - hence 
/ek ^ Cd with a proportionality factor < 1. Podolak 
et al. 1988 (hereafter PPR 88 ) argued that an upper limit 
to the proportionality factor is 1/4, though they dis¬ 
cussed the possibility for it to be smaller. Here we take 
this upper limit and assume that /ek = There¬ 

fore, the rate at which the particle gains energy due to 
frictional heating with the gas is {7r/S)CDPgR‘i 

Another source of heating is represented by the energy 
absorbed from the radiation emitted by the ambient gas 
at temperature T^, AttR^ CsCrsBTg (assuming black-body 
emission), where is the thermal emissivity of the parti¬ 
cle (here assumed a perfect black-body radiator, = 1 ^) 
and asB is the Stefan-Boltzmann constant. Similarly, 
energy is lost via radiation emitted through the particle 
surface, AirR^ CsCtsb^/- Loosely speaking, Tg represents 
the particle temperature. More precisely, as we clarify 
below, it is the temperature of an outer isothermal layer 
of the particle. 

Finally, there is energy involved in the phase transition 
of the particle’s material. If dM^/dt is the rate of change 
of the particle’s mass and all of dM^ is involved in the 
phase transition, LsdMg/dt is the energy per unit time 
absorbed (or released, depending on the sign of dMs/dt) 
in the process, where Lg is the energy per unit mass 
required to vaporize the substance. 

Accounting for all heating and cooling sources pre¬ 
sented above, the energy balance equation takes the form 


^ttRIpsCs^ = ^CopgRt |Vg - V, 


+ 47ri?2e,cTsB (T^-T^) 

dM, 


+ Ls 


dt 


( 12 ) 


In Equation (12), Cs is the specific heat of the parti¬ 
cle. Whipple (1950) estimated that, for typical mete- 
oritic material, the the left-hand side may be ignored for 
particles smaller than 0.01 cm in radius. 

In Equation (12), the variation of the particle’s inter¬ 
nal energy (the left-hand side), assumes that the temper¬ 
ature Ts is uniform throughout the volume of the body. 
Such assumption requires that there be no temperature 
gradient inside the body, i.e., that internal heat conduc¬ 
tion be infinite. This may indeed be the case for small 
particles but, as the particle radius increases, the pres¬ 
ence of a temperature gradient within the body becomes 
increasingly non-negligible. Eor example. Love & Brown¬ 
lee (1991) concluded that a significant temperature gra¬ 
dient may begin to appear across particles with a diam- 


This is a very good approximation for ice (and water), and 
typically a reasonable approximation for silicates. 


eter larger than ^ 0.1 cm when Tg ~ 1500 K. Therefore, 
the isothermality assumption advocated in Equation (12) 
may be justified only in an outer shell of the body. 

In order to evaluate the thickness of the surface layer 
of a particle, which may be approximated as isothermal 
at temperature T^, we follow the approach of McAuliffe 
& Christou (2006), based on the work of Love & Brown¬ 
lee (1991). A measure of whether or not a temperature 
gradient develops inside a heated body can be derived 
from the Biot number, which is defined as 

(13) 

where he is a characteristic heat transfer coefficient, with 
the units of an energy fiux per unit temperature, Ic is a 
characteristic length, and is the thermal conductivity 
of the material. The quantity he is intended to represent 
the rate of heat exchange between the body and the sur¬ 
rounding environment, as a function of the difference of 
temperature between them. The characteristic length le 
is typically defined as the volume-to-surface ratio. Eor a 
sphere, le is a third of the radius. 

It is customary to assume that temperature gradients 
inside a given substance are negligible for S < 0.1 (e.g., 
Lienhard & Lienhard 2008). Love & Brownlee (1991) ap¬ 
proximated the characteristic heat transfer coefficient he 
of a layer at temperature Tg as asB^f • Thus, from Equa¬ 
tion (13), one can approximate the maximum thickness 
of the isothermal layer dg to 




= 0.3 


A. 

^SbT3 


(14) 


Obviously, dg has an upper bound at Rg^ in which case 
the body can be considered as fully isothermal. 

Therefore, in general, we will assume that heating and 
cooling processes affect only a surface layer of the body, 
of thickness dg^ rather than its entire volume. In this 
approximation. Equation (12) can be re-written as 


- {Rs - 5sf]psCs^ = ^CoPgRl 


+ 47ri?2e,<7sB [T^-Tf) 


-L. 


dM, 

dt 


(15) 


Note that, in the above equation, the particle radius may 
vary with time due to ablation. At Tg = 100 K, the max¬ 
imum isothermal depth, dg, of an icy particle is few tens 
of meters, and somewhat less than ten meters at 150 K. 
Since Ag varies by a factor less than 2 between the two 
temperatures (see Table 3), this change in dg is mainly 
dictated by the increased heat exchange with the sur¬ 
roundings. The quoted depths become larger for rocky 
(quartz) bodies by a factor of ^ 3 (see Table 3), but they 
are broadly in accord with the estimate of PPR 88 , who 
concluded that heating and cooling would affect only a 
relatively thin layer of large, km-size bodies. This ap¬ 
proach, however, is rendered necessary by the fact that 
planetesimals may spend most of their time in the cool 
circumstellar and circumplanetary disk environments. 

In Equations (14) and (15), both the thermal conduc- 
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Table 3 

Material’s Properties 


Symbol 

Ice 

Rock 

Ice+Rock 

Ps^ [gcm-3] 

1.00 

2.65 

1.33 


1.00 

1.00 

1.00 

LT [ergg 1] 

2.83 X 10i° 

8.08 X lO^o 

2.83 X 10“ 

[ergg-i] 

2.50 X 10“ 

7.92 X 10^0 

2.50 X 10“ 


18.0 

60.1 

25.0 

as^ [dynecm“^] 

10® 

1- 

o 

1—1 

10® 

as [ergg-iK-i] 




at 50 K 

4.35 X 10® 

9.56 X 10® 

2.99 X 10® 

at 100 K 

8.30 X 10® 

2.67 X 10® 

6.05 X 10® 

at 200 K 

1.58 X lO’’ 

5.43 X 10® 

1.17 X 10^ 

Xs^ [ergs ^ cm ^ K 



at 50 K 

1.33 X 10® 

5.89 X 10® 

1.85 X 10® 

at 100 K 

6.41 X 10® 

2.09 X 10® 

8.33 X 10® 

at 200 K 

3.10 X 10® 

9.55 X 10® 

3.98 X 10® 


^ Density. 

^ Thermal emissivity. 

Specific vaporization energy of the solid phase. 

^ Specific vaporization energy of the liquid phase. 
^ Mean molecular weight. 

^ Nominal compressive strength at Rs = 10^ cm. 
s Specific heat. 

^ Thermal conductivity. 


3.3. Particle Ablation 

The heat deposited in the outer layer of a body can 
cause phase transitions of its material, and hence mass 
loss. Here we consider that mass loss is caused by tran¬ 
sition to the gas phase. The rate at which vaporization 
removes mass from a solid body can be approximated 
by the Hertz-Knudsen-Langmuir equation (e.g., Blottner 
1971; Campbell-Brown & Koschny 2004, and references 
therein). Indicating with Py and /i^, respectively, the 
vapor pressure and the mean molecular weight of the 
material, arguments from the kinetic theory of gases im¬ 
ply that the flux of atoms/molecules leaving the surface 
of a body is jasm}iPy/{k-BTs)Vs/^^ where Vs is the aver¬ 
age thermal speed of atoms/molecules in the vapor (e.g., 
Mihalas & Weibel Mihalas 1999): 




TT Rsmu 


(18) 


Integrating the flux over the surface of the (spherical) 
body, we have that the mass loss rate is 


dMs 

dt 


-AirRiPy 


27rkBTs 


(19) 


tivity As and the specific heat Cs are functions of We 
use piece-wise fits to the data reported by Haynes (2011) 
and Jensen et al. (1980) for ice and by Powell et al. (1966) 
and Chase (1998) for quartz (Si02). The specific energy 
of vaporization, Lg, is instead approximated as constant 
(see discussion in PPR88). Here, rather than rocks, we 
consider a mixture in which rocks are embedded in an 
icy matrix. Indicating the mass fractions of ice and rock 
with Xice and Xrock (so that Xice + Xrock = 1), the specific 
heat of the mixture is given by 

U=XiceCr + XrockCr". (16) 

The thermal conductivity of the mixture is approximated 
as 

A, = Ar [(1 - , (17) 

where 5' = XrockPT/iXroduPT + XicePU'^) is the frac- 
tion of the volume occupied by rock. The quantities 
and 4>/ represent efficiency factors for the thermal con¬ 
ductivity of a mixed medium, composed of a matrix of 
one material embedding grains of a second material (see 
discussion in Prialnik et al. 2004). Rock would consti¬ 
tute the matrix of the mixed medium for T > 0.5. In 
Equation (17), RRd are both functions of T and 
^rock/^ij^e, respectively given by Equations (25) 

and (26) of Prialnik et al. (2004). Equation (17) con¬ 
verges to for T ^ 0 (both RRd ^ 1) Rud to 

;^rock for ^ ^ 1 ^ ;^rock/;^iTe)^ 

use an ice mass fraction xice = 0.6, hence T 0.334. 
The medium remains mixed throughout the evolution 
and possible effects of differentiation (Mosqueira et al. 
2010) are ignored. A summary of some material’s prop¬ 
erties is listed in Table 3, including values of Cs and 
at three representative temperatures. 

^ Here, the material is assumed to be compact, and possible 
effects due porosity, inhomogeneity and impurity of the substance 
are neglected. 


and Py = Py{Ts). Equation (19) assumes that the vapor 
is rapidly carried away from the body’s surface, i.e., the 
partial pressure of the vapor in the gas is unimportant. 
In case of the ice-rock mixture, vapor carries away the 
icy matrix first (due to higher vapor pressure), but we 
assume that the rocky material embedded in the matrix 
is also lost by appropriately modifying 

The vapor pressure can be obtained by integrating the 
Clausius-Clapeyron equation. The resulting function de¬ 
pends on a number of constants that are fixed using phys¬ 
ical arguments and experimental data. For icy bodies, at 
temperatures below the melting point (T^ = 273.16K), 
we use the formula of Washburn (1924), which reads 

= ao P aiTs a2T^ P ^ P a^logTs. ( 20 ) 

For Py expressed in units of dyne/cm^ (= 0.1 Pa), the 
constants Oi are given in Table 4. Although Washburn’s 
formula dates back 90 years, it agrees very well with the 
2011 release of the sublimation pressure of ordinary wa¬ 
ter ice from The International Association for the Prop¬ 
erties of Water and Steam. Washburn (1924) found that 
Equation (20) satisfactorily reproduced the experimen¬ 
tal data accessible to him (above ^ 170 K). We find that 
this formula actually gives a good fit to all the values 
of the vapor pressure of ice reported by Haynes (2011), 
which extend down to 50 K. 

Above the melting point and below the critical temper¬ 
ature, Ter = 647.096 K, we use the fitting function from 
Wagner & PruB (2002) 

+ ), ( 21 ) 

where Vs = (1 — T^/Tcr), Per = 2.2064 x 10^ dyne/cm^ is 
the vapor pressure at the critical temperature, and the 
constants ai can be found in Table 4. We find that the 
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Table 4 

Constants in Vapor Pressure Formulas^ 


Equation 

do 

ai 

a2 

as 

a4 

as 

(20) 

-2445.5646 

8.2312 

-0.01677006 

0.0000120514 

-3.632266 


(21) 

-7.85951783 

1.84408259 

-11.7866497 

22.6807411 

-15.9618719 

1.80122502 

(22) 

31.82319964 

46071.4304 

58.883 





^ The pressure is in units of dyne/cm^. 


vapor pressure from Equation (21) becomes larger than 
that from Equation (20) for Tg > 272.84 K, hence we use 
this temperature value for the transition between the two 
formulas, Equations (20) and (21). 

Alternatively, the vapor pressure can be derived by fit¬ 
ting expreimental data. Eor quartz, we adopt the fitting 
function published on the Chemistry WebBook of the 
National Institute of Standards and Technology (NIST) 

lnP„ =ao-(22) 

The constants are displayed in Table 4 for Py expressed 
in units of dyne/cm^. This NIST fit applies over a limited 
range of temperatures, and it should be considered as 
an extrapolation at lower temperatures and up to the 
critical temperature. Ter = 4500 K. However, since the 
vapor pressure of the icy matrix is much higher, the mass 
loss of mixed-composition particles is also governed by 
Equations (20) and (21). 

At temperatures greater than Ter, there is no distinc¬ 
tion between the vapor and the liquid phase and the mass 
loss rate is energy limited (see PPR88). Eor T^ > Ter, 
the mass vaporization rate is 


dMs 

dt 


1 


47ri?2 (t 4 _ j.4) 




(23) 


When Ts = Ter, a particle evaporates at a constant tem¬ 
perature (Hood & Horanyi 1991). At and beyond the 
critical temperature, any net energy input is used for 
ablation. If there is a net energy output (i.e., when ra¬ 
diative cooling becomes larger than the sum of frictional 
and radiative heating) the vaporization rate is set to zero. 

Equations (19) and (23) are applied together with the 
Equations of motion (4), (5), and (6) under the hypoth¬ 
esis of isotropic mass loss, where the isotropy is with 
respect to the center of mass of the moving body. In 
other words, it is assumed that the absolute momenta 
of the escaping mass are equal to those said mass would 
have if it was attached to the moving body (Kopal 1978). 


3.4. Fracturing and Break-up of Planetesimals 

A solid body acted upon by external forces is stressed 
to some degree. In case of a spherical body, if the stress 
overcomes the compressive strength of the material, the 
body can fracture. Pollack et al. (1979) (see also Baldwin 
& Sheaffer 1971) approximated the differential force (per 
unit surface area) across a body traveling though gas as 
the dynamical pressure 

Pdy = \pg\^g-^B?- ( 24 ) 


Non-spherical bodies are also subject to bending, hence 
they can fracture at stresses lower by about an order of 
magnitude, i.e., once the tensile strength is exceeded (see 
discussion in Baldwin & Sheaffer 1971). The compressive 
strength has typically an inverse dependence on the body 
size, the body temperature, and the material porosity 
(Petrovic 2003). Material strengths are also sensitive to 
the rate of strain (e.g., Lange & Ahrens 1983), i.e., the 
rate at which the external force is applied. This may 
be especially important for large bodies. Simulations 
and laboratory experiments suggest that the strength of 
rocky, iron, and icy bodies is proportional to 1/Rs to 
a some power, which is typically between 0.3 and 0.5 
(Housen & Holsapple 1999; Benz & Asphaug 1999). 

A fractured body can quickly break apart, unless it 
is held together by its own gravity, which occurs if the 
radius exceeds 

(Pollack et al. 1979, 1986). Here we assume that if 
Rs < Rdy and the dynamical pressure exceeds the com¬ 
pressive strength of the particle’s material, the body is 
completely disrupted and the fragments quickly dissolve 
(which is probably a good approximation if the fragments 
are sufficiently small). If Rg > Rdy the body does not 
break apart, independently of Pdy 

The compressive strengths of planetesimals are largely 
unknown. Data obtained from the fragmentation of 
stony and iron meteorites in the Earth’s atmosphere im¬ 
ply strengths within the range from 10^ to 10^ dyne/cm^ 
(Ceplecha 1993; Petrovic 2001; Popova et al. 2011). The 
compressive strength of solid ice is on the order of 
lO^dyne/cm^ (Petrovic 2003), though it is expected to 
be lower for porous ice (Cox & Richter-Menge 1985). 
The inferred compressive strength of primitive icy bod¬ 
ies in the solar system, such as comets, is much smaller, 
< lO^dyne/cm^ (Toth & Lisse 2006). PPR88 argued 
that the old age of comets may have significantly altered 
their mechanical properties through outgassing. Biele 
et al. (2009) also pointed out that strengths measured 
from comets may be affected by pre-existing faulting. 
Thus, the compressive strength of “young” icy plan¬ 
etesimals may as well be in the range from r\j 10® to 
~ 10® dyne/cm^. 

Given the large uncertainties, we set the material com¬ 
pressive strength to ag^/lkm/Rg (Holsapple 2009) . The 
nominal strength, (jg, at the Ikm-scale size is 10^ and 
lO^dyne/cm^ for icy and rocky planetesimals, respec¬ 
tively (see Table 3). This choice of the compressive 
strengths implies that, according Equation (25), only 
icy (rocky) bodies whose radius is smaller than 10 km 
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Table 5 

Grid Structure 


Level 

Nr 

Ne 



Volume^ 


1 

243 

22 

423 



disk 



2 

84 

24 

84 

9.13 

X 

2.53 

X 

9.15 

3 

104 

34 

104 

5.65 

X 

1.79 

X 

5.66 

4 

124 

44 

124 

3.37 

X 

1.16 

X 

3.37 

5 

164 

64 

164 

2.23 

X 

0.84 

X 

2.23 

6 

244 

84 

244 

1.66 

X 

0.55 

X 

1.66 

7 

404 

104 

404 

1.37 

X 

0.34 

X 

1.37 

8 

724 

144 

724 

1.23 

X 

0.24 

X 

1.23 


^ Volume is in units of except for grid level 1. 


(« 20 km) can fragment, if the dynamical pressure is just 
marginally larger than the compressive strength. Larger 
bodies can fracture more easily, but are held together 
by their own gravity. However, the maximum radius for 
break-up increases as Pdy increases. 

4. NUMERICAL METHODS 
4.1. Solution for the Disk Hydrodynamics 

The disk’s gas is described as a continuum viscous 
fluid via the Navier-Stokes equations (see, e.g., Mihalas 
& Weibel Mihalas 1999), written in terms of the spe¬ 
cific linear momentum and the specific total angular mo¬ 
menta of the gas in a rotating frame (see D’Angelo et al. 
2005). These equations are solved in a stepwise fashion 
by means of a finite-difference code. The solution of the 
advection term, referred to as the transport step, applies 
the monotonic transport of van Leer (1977) and uses an 
operator-splitting technique (see Stone & Norman 1992) 
to cope with the three dimensions. In the source step, 
the other terms of the equations are taken into account, 
namely the apparent forces, the gradients of pressure and 
gravity, and the viscous stresses. Numerical stability is 
ensured by constraining the integration time step. At, 
according to the Courant-Friedrichs-Lewy condition (see 
Stone & Norman 1992). Overall, the algorithm is second- 
order accurate in space and effectively second-order ac¬ 
curate in time (e.g.. Boss & Myhill 1992). The code was 
compared against other fluid dynamics codes in studies 
involving problems of tidal interactions between planets 
and disks (de Val-Borro et al. 2006; Masset et al. 2006; 
de Val-Borro et al. 2007). 

The Navier-Stokes momentum equations are dis¬ 
cretized over a spherical polar grid with constant spac¬ 
ing in all three coordinate directions. The code allows 
for grid refinements by means of a nested-grid technique 
(D’Angelo et al. 2002, 2003a). The increase of volume 
resolution is a factor of 2^ for any level added to the 
grid system. In this study, we employ a grid system with 
8 levels, the details of which are given in Table 5. The 
first level encloses the entire disk, whereas additional lev¬ 
els enclose smaller and smaller disk portions around the 
planet. In Table 5, Nq^ and indicate the number 
of grid points along the correspondent coordinate direc¬ 
tions. The last column gives the volume occupied by each 
grid level, where the lengths are in units of Rh- Overall, 
the grid system contains about 103 million grid elements. 
The region of the wider circumplanetary disk, typically 
taken as ^ Rh/ 4 around the planet, is discretized over 
the 8^^ grid with more than 12 million gird elements. 


The spatial resolution on the first grid level is such that 
Arjop Op A(l)/{rsm0) 0.015 and apAOjr 0.013. 

On the 8^^ grid level, the linear resolution around the 
planet is ~ 10' ^Up ^ 1.4x 10 ^ Rh, which is about equal 
to Jupiter’s current radius, Rj, at 5.2 AU. Note that the 
actual radius of the planet, Rp, at these late stages of 
accretion (i.e., when it is no longer accreting substantial 
quantities of gas compared to its mass) is likely > 1.3 Rj 
and < 1.8 Rj (Lissauer et al. 2009). We adopt the value 
Rp = 1.6 Rj. 

We apply boundary conditions at the inner and outer 
disk radii, rmn and rmx, using the procedures of de Val- 
Borro et al. (2006). Boundaries at the disk surface {0 = 
^mn) aud at the equatorial plane are handled as in Masset 
et al. (2006). In these calculations, we do not account for 
accretion on the central star which, in conjunction with 
accretion on the planet, can alter the density in the disk 
interior of the planet’s orbit (Lubow & D’Angelo 2006). 

4.2. Solution for the Planetesimal Thermodynamics 

The set of Equations (4)-(6) is completed by the equa¬ 
tions to obtain the spherical polar coordinates of a par¬ 
ticle 

dr 


_ L(j) _ Q 

dt {rsinOy 

The system of first order ordinary differential equations 
(ODE), represented by Equations (26), (4), (5), (6), (15), 
and (19) or (23), is solved numerically by means of a vari¬ 
able (arbitrarily high) order and variable step-size Gragg- 
Bulirsch-Stoer extrapolation algorithm (Hairer et al. 
1993). Indicating with At the time step of the hydro- 
dynamical calculation at time t (see Section 4.2), the 
system of ODE is integrated three times, according to 
the step-size sequence (At/4, At/2, At/4), using gas field 
distributions centered at (t, t+At/2, t-FAt), respectively. 

The ODE solver chooses automatically the order of the 
algorithm and a series of appropriate internal time inter¬ 
vals so to advance the solution to the required end time. 
The algorithm’s order and the length of each time in¬ 
terval are constrained by user-supplied tolerances on the 
local truncation error of the solution, which is estimated 
from the comparison of solutions at different orders. Here 
we apply tolerances in the range from 2 x 10“^^ to 10“^^. 

The 3D gas field distributions of T^, and v^, are 
interpolated in space at the position of the particles by 
means of a second-order accurate algorithm, based on 
monotonic harmonic means (van Leer 1977), according 
to the approach of D’Angelo et al. (2002), extended to 
three dimensions. The advantage of this method rests on 
its capability of handling discontinuities and shock-like 
conditions in the gas. The spatial as well as temporal in¬ 
terpolations are performed on the gas field distributions 
with the highest available resolution, which are those cal¬ 
culated on the most refined grid level where the particle 
is located. 

Several tests of the planetesimal thermodynamics 
solver are presented in Appendix B. These include stan¬ 
dard two- and three-body problems, drag-induced orbital 
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decay and free-fall of particles, and various thermal evo¬ 
lution problems. 

4.3. Gas and Particle Accretion 

The accretion of gas onto a gap-opening planet is a 
complex problem. It was suggested by D’Angelo et al. 
(2003b) and Bate et al. (2003), and later confirmed (see 
the recent studies by Tanigawa et al. 2012; Ayliffe & Bate 
2012; Gressel et al. 2013; Szulagyi et al. 2014, and refer¬ 
ences therein), that gas mostly proceeds off the mid-plane 
(at and above the surface) of the disk around the planet, 
prior to accreting on its envelope. We do not model the 
planet’s envelope here and adopt a prescription for gas 
accretion along the lines of D’Angelo et al. (2003b) and 
D’Angelo & Lubow (2008), a procedure that is direction- 
ally unbiased. The spherical volume around the planet 
from which gas is removed to mimic accretion extends for 
lARp in radius (~ 2.2i?j), which makes the procedure 
independent of the mode of gas delivery to the planet’s 
envelope. At such short distances, the thermal energy 
of the gas is much smaller than the gravitational energy 
binding the gas to the planet (Bodenheimer & Pollack 
1986), hence gas cannot escape. Planet formation calcu¬ 
lations do indicate that the rates of accretion calculated 
in this manner correspond to the actual rates of envelope 
growth (Lissauer et al. 2009). 

The accretion of planetesimals is a simpler problem, 
since arbitrarily close encounters with the planet are al¬ 
lowed. During close approaches, the gravitational poten¬ 
tial of the planet is always used, and no regularization is 
applied (see, e.g., Bodenheimer et al. 2006). We adopt 
two criteria for accretion: if the particle approaches the 
planet within a distance < Rp (1.6 i?j), a head-on impact 
is assumed; otherwise, and if the distance of approach is 
< 2.2 Rp (3.5 i?j), the particle is deemed as accreted if 
its relative velocity is less than the escape velocity from 
the planet at that distance. 

5. DISK DENSITY AND DYNAMICS 

D’Angelo & Marzari (2012) performed calculations of 
circumstellar disk evolution driven by viscous diffusion 
and photoevaporation, exploring ranges of stellar EUV 
luminosity, initial disk mass, initial mass distribution, 
and gas kinematic viscosity representative of the proto¬ 
sun and the early solar nebula. These parameters can 
be constrained by the requirements that the disk’s gas 
lifetime be shorter than 20Myr (e.g., Haisch et al. 2001; 
Pascucci et al. 2006; Roberge & Kamp 2011; Williams & 
Cieza 2011; Bell et al. 2013) and longer than the forma¬ 
tion time of a Jupiter-mass planet at ~ 5 AU, which is 
> 1 Myr (Hubickyj et al. 2005; Alibert et al. 2005; Lis¬ 
sauer et al. 2009; Movshovitz et al. 2010; Mordasini et al. 
2011). 

In Eigure 1, we plot results from some models of 
D’Angelo & Marzari (2012), in particular the surface 
density at 5 AU versus time. The ratio of the initial disk 
mass to the stellar mass is indicated in the top-right cor¬ 
ner. The curves indicate that, after ^ 1 Myr, the density 
is at most ~ 200gcm“^, and is typically smaller than 
100gcm“^ (but could be much smaller, « 20gcm“^). 
After ^ 2 Myr, the surface density ranges from ^ 10 to 
^ 50gcm“^. In all the models in the figure, the initial 
surface density inside of about 10 AU is proportional to 



Time [Myr] 

Figure 1. Surface density at 5AU versus time in disk models 
whose evolution is driven by viscous diffusion and photoevapora¬ 
tion. The numbers in the legend indicate the initial disk mass, in 
units of M*, within 40 AU of the star. The models also employ 
different gas kinematic viscosities and EUV fluxes emitted by the 
star. See text for an explanation of the shaded area. Data from 
the models of D’Angelo &: Marzari (2012). 

l/yY, consistent with that used in the hydro dynamical 
calculations. 

Therefore, given the range of possible disk densities at 
the time of Jupiter’s formation (between, say, ^ 1 and 
^ 4Myr), we consider two values for the unperturbed 
surface density, Eq, at 5.2 AU of ~ 10 and ~ 100gcm“^, 
which correspond to mass densities po ~ 10“^^ and po ~ 
10“^^gcm“^. The shaded area in Eigure 1 represents 
the area covered by our choice of Eq. The disk mass in 
units of M^, inside of « 21 AU, is ^ 8 x 10 “^ Eq, where 
Eo is in units of gcm“^. 

Tidal interactions between the disk’s gas and the 
planet excite density waves at Lindblad resonances (Gol- 
dreich & Tremaine 1980) and deplete the gas within a few 
Rh from the planet’s orbit, where tidal torques exceed 
viscous torques (Lin & Papaloizou 1986). Residual gas 
is still present in the tidal gap region (as shown below), 
even at a viscosity much smaller than that adopted here. 

The main features of the surface and volume density, 
on a global disk scale (> a^), are illustrated in Eigure 2. 
In all cases, densities are normalized to either Eq or po- 
Both spiral density waves and the tidally-produced gap 
are visible in the top-left panel, while the plot on the 
right shows, more quantitatively, the volume density in 
the disk’s equatorial planet, at several distances from the 
planet’s orbit. The residual gas in the tidal gap region 
is also visible as a function of the azimuthal angle. The 
center panels illustrate the density on orthogonal disk 
slices passing through the planet’s position, on different 
length scales. The surface density of the region around 
the planet’s Roche lobe is shown in the bottom panels 
(the left panel also shows the Roche lobe trace and the 
positions of the Lagrange points Li - lower cross - and 
1/2 - upper cross). Radial cuts of E^/Eq at various az¬ 
imuthal angles are plotted in the right panel. All images 
are saturated in order to improve the contrast between 
low and high density regions. 

The rotation curve of the unperturbed disk is affected 
by the gas pressure gradient, which depends on both den¬ 
sity and temperature gradients. In terms of the Keple- 
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Figure 2. Top-left. Color scale rendering of log(S^/So) as a function of the azimuthal angle from the planet and the distance from 
the star in units of ap. Top-right. Volume density, pg in the disk’s mid-plane, normalized to po as a function of the azimuth at various 
distances from the planet’s orbit, as indicated in the legend. Center. Color scale rendering of log(p^/po) on disk slices passing through 
the planet. The latitude z/r) is on the vertical axis and the linear distance from the star (left and middle) or angular distance from 
the planet (right) is on the horizontal axis. Bottom. Color scale rendering of log (S^/Sq) in the proximity of the planet’s Roche lobe (left) 
and S^/So versus radial distance (right), at various separation angles, \4> — 0p|/7r, as indicated. Thicker lines are for > 0^. 


rian velocity, tk, the (absolute) azimuthal velocity of the 
gas, in absence of the planet and in the mid-plane of the 
disk, would be (e.g., Tanaka et ah 2002) 



UK 



(27) 


which accounts for the fact that the unperturbed density 


is pg oc 1/r^/^ and Tg (x 1/r. The subscript “tx” stresses 
the fact that this expression does not account for the 
perturbation induced by the planet. Velocity )u dif¬ 
fers by less than 1% from the Keplerian velocity. The 
ratio of the perturbed velocity, , to that in Equa¬ 
tion (27) is shown in Figure 3. The largest deviations 
from the unperturbed rotation curve occur around the 
edges of the gap, where the magnitude of the density 
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Figure 3. Absolute azimuthal velocity of the gas, averaged over 
27r in azimuth around the star and normalized to the unperturbed 
velocity of Equation (27). The radial pressure gradient of the gas 
induces sub/super-Keplerian rotation at the inner/outer gap edge. 


gradient is the largest (see top-left panel in Figure 2). 
The negative/positive pressure (radial) gradient triggers 
a sub/super-Keplerian rotation around the inner/outer 
edge of the gap. In fact, centrifugal balance requires 
that = '^K + !Pg)^Pg !using Equation (2), 


V 


A 

0 


'^K 



r J \ d\nr 



(28) 


which reduces to Equation (27) for the unperturbed disk 
case. Eigure 3 and Equation (28) suggest that the magni¬ 
tude of the gradient dlnpg/dlnr is marginally larger at 
the outer edge of the gap than it is at the inner edge. If a 
particle moved at a Keplerian speed, on average it would 
experience a tail wind when orbiting near the outer edge 
of the density gap and a head wind when orbiting near 
the inner gap edge. 


5.1. Circumplanetary Disk Thermodynamics 

Le us introduce a local reference frame {0';x^y^z}, 
with origin O' on the planet, coordinate x pointing away 
from the star, y pointing toward the direction of orbital 
motion, and 2 ) pointing away from the disk’s equatorial 
plane {0 = tt/2) in the direction 0 = 0. 

The normalized density, pgjpo^ is illustrated in Eig¬ 
ure 4, on length scales < on vertical (top) and equa¬ 
torial (bottom) slices passing though the position of the 
planet. The density contour levels in the lower panels 
of the figure indicate that a perturbation in the form 
of a spiral wave propagates toward the planet, but it 
does not propagate in the inner disk, within a distance 
of 0.1 or about 75 Rj of the planet. This is in agree¬ 
ment with previous 3D calculations. Note that while all 
of Jupiter’s regular satellites orbit inside 30 Rj of the 
planet, irregular satellites have semi-major axes between 
^ 100 -Rj ~ 0.13-Rh nnd ^ 425 Rj ~ 0.56 Rh (c.g., Je- 
witt et al. 2004; Jewitt & Haghighipour 2007). Beyond 
^ 240 Rj ~ 0.32 Rh, these satellites are all on retrograde 
orbits. 

Line plots of the vertical density distribution are shown 
in the top panel of Eigure 5, at various distances from 
the planet’s location (see figure caption for details). The 
volume and surface densities within the inner Roche lobe 
are shown in the middle and bottom panels, respectively. 
A close-up of these quantities around the planet can be 
seen in the insets. The surface density between 0.02 
and 0.15 Rh of the planet roughly declines as a —0.3 


power of the distance. This slope becomes approximately 
— 1 between ^ 0.15 and 0.5 Rh- 

The two-dimensional models of D’Angelo et al. (2003a) 
showed that circumplanetary disks are typically thick, 
with an aspect ratio ranging from ^ 0.2 to ^ 0.4, depend¬ 
ing on the thermal state of the disk. This conclusion was 
later confirmed by several other studies (e.g., Machida 
et al. 2008; Ayliffe & Bate 2009; Martin & Lubow 2011; 
Tanigawa et al. 2012). The top-right panel of Eigure 4 
also predicts a thick disk, in agreement with previous 
results. A direct measurement from the data in the fig¬ 
ure, of the location where there is a sharp density drop, 
yields an aspect ratio of about 0.36 within « 0.13 Rh of 
the planet. 

Although we use a local isothermal equation of state 
(Equation (2)), which becomes effectively isothermal in 
the small region around the planet occupied by the cir¬ 
cumplanetary disk, we set a gas temperature using simple 
arguments based on local heating via viscous dissipation 
and black-body radiation from ambient gas, and verti¬ 
cal radiative cooling. Indicating with -h the 

effective disk temperature, Te, is given by (Pringle 1981) 


rp4: rpA 

^ e -^n 


3 GMpMp 
Stt (Jsb 



(29) 


where refers to the cricumstellar disk tempera¬ 
ture in Equation (3) and Rp = 1.6 Rj. The max¬ 
imum of the right-hand side of Equation (29) is ~ 
0.00677 GMpMp/( cfsbRI) and occurs at f = (49/36) Rp. 
The gas accretion rate on the planet, measured from the 
calculation (see Section 4.3), is Mp ^ 2 x 
However, the accretion rate involved in Equation (29) 
is actually that through the disk, which is only a small 
fraction of Mp (as mentioned in Section 4.3). Although 
this fraction varies with distance from the planet, based 
on the analysis of Tanigawa et al. (2012), we simply ap¬ 
proximate it to about 1/6. Since there is no physical 
boundary at R^, the interface between the disk and the 
planet, for f^ < R^ we set Tg equal to 600 K, the 

effective temperature of a protojupiter right after most 
of the envelope has been accreted^ (Lissauer et al. 2009). 
We neglect possible heating effects in the circumplane¬ 
tary due to irradiation by the planet. 

Eollowing Lunine & Stevenson (1982), the vertical tem¬ 
perature can be derived by considering energy transfer 
via radiation in the vertical direction, which leads to the 
the equation 

^ = -3p,k^{t:-T^) (30) 

where /dr is a frequency-integrated opacity, which we 
approximate as the Rosseland mean opacity, and H ^ 
0.36 r is the local circumplanetary disk scale height (see 
top panels of Eigure 4). By solving Equation (30), we 


® In the models of Lissauer et al. (2009), the effective temper¬ 
ature of a protojupiter, after most of the gaseous envelope has 
been acquired, may initially depend on the accretion history of the 
planet, but it subsequently converges to ~ 600 K within a few 10^ 
years (see their Figure 10) and to ~ 500 K by lOMyr (Marley et al. 
2007). 
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Figure 4. Color scale rendering of log (p^/po) on disk slices passing through the planet’s location. Vertical slices are illustrated in the top 
panels whereas slices in the equatorial plane are illustrated in the bottom panels. All distances from the planet are in units of the planet’s 
Hill radius, Rn ~ 0.0689 a,p. 


have 

R = R + {R - R) j , ( 31 ) 

with a height-dependent optical depth tr = pgKYi{H — 
\z\), assumed to be non-negative. The quantity pg is a 
vertically averaged volume density between heights |^| 
and H. Although the equation above is derived for a 
vertically constant opacity, we assume that /^r is either 
constant or a linear function of Tg. Thus, the solution of 
Equation (31) typically requires a root-finding iteration 
procedure, for which we use an algorithm based on the 
Brent’s method (Brent 1973). 

We recall that, in Equation (31), Tg indicates the tem¬ 
perature in the circumplanetary disk, is the tempera¬ 
ture in the circumstellar disk (Equation (3)), and T^—T^ 
is given by Equation (29). At | 2 (| = we have that 
Tg = Tg whereas, for |^| > i7, we impose an exponential 
decline over height of Tg until it matches Tn. 

The underlying local isothermal approximation used in 
the hydrodynamics calculations implies that the opacity 


of the medium (gas+dust) is low. If we apply a gas- 
dominated opacity, /^r = 10 “^cm^g“^, as contemplated 
in some of the models of Canup & Ward (2002), then Tg 
becomes about equal to the effective temperature Tg, as 
can be seen by comparing the red and blue curves in the 
two panels of Figure 6. 

In Equation (31), as anticipated above, we use an opac¬ 
ity that is a piece-wise function of the temperature. Be¬ 
low the vaporization temperature of dust grains (Pol¬ 
lack et al. 1994), we set /dr = O.Olcm^g”^, which may 
be the case in an evolved disk where grains have under¬ 
gone significant growth (e.g., D’Alessio et al. 2001). We 
assume that /dr linearly transitions to a gas-dominated 
opacity, /dr = 10 “^cm^g“^, in the temperature interval 
from 1600 K to about 2000 K, and it remains constant at 
larger temperatures. For this opacity law, the equatorial 
temperature of the circumplanetary disk is illustrated as 
a black line in Figure 6, for both reference background 
densities po = 10“^^ (top) and 10“^^gcm“^ (bottom). 
At distances from the planet f > 0.2 the radial dis¬ 
tribution of temperature merges with the temperature 
distribution in the circumstellar disk. 

Since Equation (31) is not consistent with the equa- 
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Figure 5. The top panel shows the normalized volume density, 
pgjpo^ along the vertical direction at various distances from the 
planet. From profiles with lower to higher peak density, the dis¬ 
tance is 0.5, 0.35, 0.2, 0.1, 0.05, and 0.02 The middle and 
bottom panels show, respectively, the normalized the (mid-plane) 
volume and surface density around the planet. The top a^-axis of 
each panel gives the distance in units of the planet’s semi-major 
axis, Up, while the bottom axis is in units of Rn. 


tion of state applied in the hydrodynamics calcula¬ 
tions (Equation (2)), we also consider cases in which 
the gas temperature is given everywhere by Equa¬ 
tion (3). In these calculations, the circumplanetary disk 
is basically isothermal with a gas temperature 
{Rgmi{/kB){GM^/ap){H/ap)‘^. The lower temperature 
close to the planet may affect the ablation history, and 
hence the mass evolution, of some planetesimals. 



Figure 6. Temperature in the equatorial plane of the circum¬ 
planetary disk, from Equation (31), as a function of the distance 
from the planet, for background volume density po = 10“^^ gcm“^ 
(top) and 10“^^gcm“^ (bottom). The shaded areas indicate 
distances r < Rp. The black curves use a Rosseland mean 
opacity = 10“^cm^g“^ for Tg < 1600K, which decreases 
linearly toward a gas-dominated opacity (/^r = 10“"^ cm^ g“^). 
The red curves are for a constant, gas-dominated opacity rr = 
lO^'^cm^g”^, while the blue curves represent the effective disk 
temperature, Te (Equation (29)). 


6 . EVOLUTION OE PLANETESIMALS IN THE 
CIRCUMSTELLAR DISK 

We follow the evolution of planetesimals initially 
equally distributed in four size bins with radii = O.I, 
I, 10, and 100km. The planetesimals are placed on ellip¬ 
tical orbits about the star. The initial orbital eccentric¬ 
ity, eg, and inclination, U? aro randomly selected within 
the range from 0 to 0.05 and from 0 to 0.05 radiants 
2.9°). The initial argument of periapsis, longitude of 
the ascending node, and true anomaly are chosen ran¬ 
domly between 0 and 27r. At the beginning, three re¬ 
gions in semi-major axis, a^, are populated: from 0.77 
to 0.82 Up, from 0.965 a^, to 1.035 a^, and from 1.2 to 
1.25 Up. These regions are all inside 4 Rh of the planet’s 
orbit, the classical “feeding zone” for accretion of solids 
(e.g., Greenzweig & Lissauer 1990; Lissauer & Stewart 
1993). Planetesimals deployed in the corotation region 
of the planet, ap^RH/ 2 , have near-circular obits and are 
deployed with an azimuth such that |0 — 0p| > 7 r/ 3 , i.e., 
in between the triangular Lagrange points L 4 (leading) 
and L 5 (trailing). These three regions are each popu¬ 
lated with 132000 planetesimals (for a total of 0.65 or 
0.86 Mars masses of solids, depending on the material). 
Initial surface densities of solids are between ^0.05 and 
~ 0.12 gcm“^, depending on the region and the ma¬ 
terial. It is important to stress that, given the equal 
number densities per size bin, the solid mass is almost 
entirely carried by the largest bodies. Icy and mixed- 
composition planetesimals are considered in separate cal- 
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culations. The results presented here can be rescaled by 
the initial surface density of solids, provided that inter¬ 
actions among planetesimals can be neglected (see Sec¬ 
tion 6 . 1 ). 

In a planet-less disk, the orbital evolution of planetes¬ 
imals would be dictated only by gas drag (and stellar 
gravity). Thus, orbital eccentricity and inclination would 
be damped on a timescale Tdamp ^ Ws — ^g\/\^D\- Ir R 
nearly Keplerian disk, with negligible radial velocity and 
with azimuthal velocity given by Equation (27), the ap¬ 
proximation I Vs — v^l ^ (Se^/S + 2^/2 + ^^/4)^/^asflK 
can be adopted, where = (5/2) (i7/r)^ and fix is the 
Keplerian orbital frequency of the planetesimal about 
the star (e.g., Adachi et al. 1976; Ogihara & Ida 2009). 
Therefore, 


1 


^KTdrag 




y + + (32) 


The timescale for the removal of orbital energy, and 
hence for the variation of the planetesimal semi-major 
axis, is much longer, of order Tdrag/-?^ (see also Equa¬ 
tion (B5)). Our initial conditions {cs and is < 0.05) 
would lead to damping timescales, at ^ 5 AU and for 
po = 10“^^gcm“^, Tdrag > 20 orbits for r\j 0.1 km-size 
bodies {Cd ~ 6 ) and ^ 2 x 10 ^ orbits for ^ 100 km- 
size bodies {Cd ^ 0.4). It is worth noticing, however, 
that Equation (32) assumes that pg is constant along 
the trajectory of the planetesimal. The timescales for 
drag-induced orbital decay would be over two orders of 
magnitude as long. The planetesimal evolution presented 
here lasts for 580 planet’s orbits, or « 7000 years. Al¬ 
though the initial orbits of planetesimals are arbitrary, 
by the end of the calculations the spatial distributions 
of the solids are in a state of quasi-equilibrium, in the 
sense that they vary slowly over tens of orbital periods 
of the planet. Over much longer timescales, the lack of a 
supply of planetesimals from other regions (besides those 
considered here) likely inhibits a state of true dynamical 
equilibrium. 

The distributions of some instantaneous (osculating) 
orbital elements are presented in Figures 7 and 8 , at the 
end of the evolution (see the figure captions for further 
details). The histograms refer to bodies deployed in the 
regions interior and exterior of the planet’s orbit (Fig¬ 
ure 7) and in the planet’s corotation region (Figure 8 ). 
Bodies that move beyond the disk boundaries (see Sec¬ 
tion 2 ) are removed from the calculations and excluded 
from the analysis. Between « 11% and ~ 14% of the 
initial mass in solids moves out of the computational do¬ 
main by the end of the calculations, where the larger 
percentage is generated by cases with po = 10 “^^ gcm“^ 
(icy and mixed-composition bodies produce similar frac¬ 
tions). The average rates of mass loss through the bound¬ 
aries, between 2% and 3% of the initial mass per 100 or¬ 
bits, are consistent with those measured over the last 100 
orbits of the calculations. The probability of ejecting out 
of boundaries bodies of ^ 10-100 km in radius is similar 
(within factors of order unity), and higher or somewhat 
higher (depending on po) than the probability of ejecting 
< 1 km bodies. 

During the evolution, only about 0.3% of the initial 
mass is ablated. The rates of mass lost to ablation 
are roughly steady during the course of the calculations. 


Higher gas densities produce slightly more ablation than 
do lower densities. Models that use everywhere the 
gas temperature given by Equation (3), instead of the 
locally modified temperature discussed in Section 5.1, 
yield similar fractions for the ablated mass (see Sec¬ 
tion 7.1). The small percentage of ablated mass masks 
the fact that the total mass in solids is basically car¬ 
ried by 100 km bodies, which shed relatively little mass. 
Smaller bodies, however, are more prone to ablation. In 
fact, Rs ~ 10 km bodies lose ^ 5% of their total ini¬ 
tial mass and Rg ^ 1 km planetesimals shed ^ 20 % of 
their total initial mass. The fraction becomes ^ 40% for 
Rs ~ 0.1km bodies. This outcome can be understood 
from Equation (19), assuming bodies of equal tempera¬ 
ture (and hence vapor pressure, P^), which yields a ratio 
of the ablated to the total mass proportional to 1/Rs- 
Had all size bins contained equal masses, the swarm 
would have lost to ablation over 10 % of its original mass. 

Negligible fractions of the initial mass are lost through 
fracturing and break-up: ^ 10 “^ for po = 10 “^^gcm“^ 
and ^ 10“^ for po = 10“^^ gcm“^. Bodies can break up 
only if their radius is smaller than Pdy (Equation (25)) 
and the body’s compressive strength is exceeded by the 
dynamical pressure, -Pdy, which is proportional to the 
gas density (Equation (24)). Break-up of planetesimals 
occurs within 0.2 Ph of the planet, when they impact 
with the dense regions of the circumplanetary disk. 

Comparisons of the histograms show only marginal dif¬ 
ferences between the orbital elements of icy and mixed- 
composition bodies, for the same value of the gas density 
Po. As expected from the discussion above, calculations 
with Po = 10 “^^ and 10 “^^ gcm“^ provide similar distri¬ 
butions of semi-major axes, and distributions of eccen¬ 
tricities and inclinations differing mainly toward large 
values {eg ^ 0.5 and ig ^ 5°). Regardless of the plan¬ 
etesimal composition and background gas density, by the 
end of the calculations, about 5% of the remaining mass 
of solids initially placed interior of the planet’s orbit is 
scattered outside the orbit. Toward the end of the cal¬ 
culations, the average scattering rate is about 1 % of the 
remaining mass per 100 orbital periods. Around 11 % 
of the available mass of solids initially present exterior 
of the planet’s orbit is scattered inside, with an average 
scattering rate (toward the end) of about 6 % of the re¬ 
maining mass per 100 orbits. This would amount to an 
average of ^ 10 ^^ g of solids scattered toward the inner 
disk during a planet’s orbital period, if the initial surface 
density of planetesimals between 1.2 and 1.25 was 
lgcm“^. Although scattering involves planetesimals of 
all sizes, larger size objects {Rg > 1 km) are typically 
scattered more efficiently, in either radial direction, than 
are smaller size objects. Moreover, Rg ^ 0.1km bod¬ 
ies are more easily scattered outward, from inner disk 
regions, at the lower (rather than at the higher) value 
of the reference density, po. While not visible in the 
histograms of Figure 7, because of the large bin size, dis¬ 
tributions with finer sampling in semi-major axis show 
several dips, in proximity to the position of the 3:2, 5:3, 
and 2:1 mean-motion resonances with the planet, and 
to the corresponding resonant locations exterior of the 
planet’s orbit. 

Most of the mass deployed in the corotation region 
(up ^ Rh/ 2 ) remains within ^ Rn throughout the cal- 
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Figure 7. Distributions of semi-major axis (left), eccentricity (center), and inclination (right) for po = 10 (upper pair of rows) and 
10“iigcm“^ (lower pair of rows). Odd (even) rows illustrate the distributions of planetesimals deployed in the region inside (outside) 
of the planet’s orbit. The cross-hatched histograms show the initial distributions of the mixed-composition bodies. The color-shaded and 
line histograms represent the distributions, respectively, of the mixed-composition population and of both the icy and mixed-composition 
populations, about 580 orbits after deployment. 


culations (see Figure 8 ). Only a small fraction, ^ 0.5%, 
of the initial mass is scattered toward the inner disk and 
a fraction of a few percent is scattered outward. Taking 
into account all possible fates for the solids, the mass 
depletion rate of the radial region ^ Ru is around 
2 % of the initial mass per hundred orbits of the planet. 
Note that gas drag in the region along the planet’s or¬ 
bit is reduced (except very close to the planet) due to 
the density gap (see Figure 2 ). Within factors of order 
unity, the ejection probability toward the inner disk and 
that toward the outer disk are constant in the size range 
0 . 1 km ^ Rs ^ 100 km. Along the planet’s orbit, the 
largest number densities in the frame corotating with the 
planet occur around the L 4 and L 5 Lagrange points. Al¬ 
though the longitude relative to the planet of these points 
can be affected by gas drag, they effectively lie 60° ahead 
(L4) and behind (L5) the planet (for Rg ^0.1 km) owing 


to the presence of the density gap (Peale 1993). In the 
radial region of tadpole orbits, ^ 0.74i?H (e.g., Mur¬ 
ray & Dermott 2000), the number densities within a 15° 
longitude of either point are roughly constant for plan¬ 
etesimals of all sizes. The smallest number density (in 
the corotating frame) is around the collinear I /3 point. 

The gap in the planetesimal disk, due to the torques ex¬ 
erted by the planet on the solids, is illustrated in Figure 9 
as a function of the semi-major axis, for various planetes¬ 
imal radii (see figure caption for further details). The 
distributions comprise icy and mixed-composition bod¬ 
ies, initially deployed interior and exterior of the planet’s 
orbit. The positions of the gap edges move further away 
from the planet’s orbit as Rg reduces because of gas drag 
effects (see below). Differences with respect to the gas 
density appear marginal for radii Rg ^ 1 km, but they 
become more pronounced at radii Rg < 0 . 1 km (parti- 
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Figure 8. As in Figure 7, but for bodies initially placed in the planet corotation region. Top (bottom) panels refer to po = 10 (10 

gcm“^. 
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Figure 9. Semi-major axis distributions of planetesimals in prox¬ 
imity of the planet’s orbit. The histograms include bodies ini¬ 
tially placed inside and outside of the orbit and of both com¬ 
positions. The reference gas density is po = 10“^^ (top) and 
10“iigcm“^ (bottom). Histograms of different colors progres¬ 
sively include larger and larger bodies: Rs < 0.1 (lowest-count 
histogram), 1, 10, and 100km (highest-count histogram). In the 
bottom panel, Rg < 0.1km bodies do not appear because they are 
located farther from the planet’s orbit. 


cles of these sizes do not appear in the bottom panel of 
Figure 9). For po = the outer gap edge 

of Rs ^ 0.1km bodies recedes at r ~ 1.25 a^, where 
the (azimuthally) averaged rotation velocity of the gas 
exceeds the azimuthal velocity of particles because the 


radial pressure gradient of the gas is locally positive (al¬ 
though only marginally, this effect is already present at 
r ^ 1.6 Up). The inner gap edge of 0.1 km bodies is found 
at r ^ 0.67Up. To address the behavior of small frag¬ 
ments, in Section 8 we discuss some experiments con¬ 
ducted with cm-to-m size particles. 

In fact, the super-Keplerian rotation at the outer edge 
of the density gap, mentioned in Section 5 (see Fig¬ 
ure 3), can lead to planetesimal segregation, by halt¬ 
ing or pushing outward bodies that reach those radial 
locations (see, e.g, Ayliffe et al. 2012, and references 
therein). Neglecting collisions and gravitational encoun¬ 
ters among bodies, which would redistribute bodies and 
likely work against segregation, the efficiency of this pro¬ 
cess depends on the competing effects of the gravitational 
torque (perpendicular to the orbital plane) exerted by the 
planet and the vertical component of the gas drag torque 
VsXFd- Considering planetesimals on near-Keplerian or¬ 
bits, rsX(v^ —Vg) has vertical component ^ as{v^ — vk)^ 

where is given by Equation (28) (see also Figure 3), 
which is positive/negative at the outer/inner gap edge. 
For large enough objects {Cd ^ 1), the gas drag torque 
is then oc R^^pgasV^ oc R^pg whereas the gravitational 
torque due to the planet is proportional to the body 
mass and hence to yielding a segregation efficiency 
oc pgjRg. Therefore, for a given value of the local den¬ 
sity (and distance from the planet’s orbit), one can ex¬ 
pect segregation of smaller planetesimals to be more effi¬ 
cient than that of larger bodies. As expected, we do not 
observe strict segregation of planetesimals. However, it 
does appear that Rg ^ 0.1km bodies, initially placed in 
the region exterior of the planet’s orbit, are less likely to 
move toward the inner disk than are Rg ^ 1 km bodies, 
by a factor of ^ 10 (see Figure 9, top panel). Addition¬ 
ally, Rg ^0.1 km bodies are less prone to cross from the 
outer to the inner disk at the higher value of po than they 
are at the lower density, again by a factor of order 10. 

Figure 10 shows two-dimensional histograms of eccen¬ 
tricity, inclination, size, and temperature of planetesi- 
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Figure 10. Distribution of icy and mixed-composition particles 
versus semi-major axis and versus orbital eccentricity, orbital incli¬ 
nation, radius, and temperature, as indicated on the vertical axes. 
Bodies initially placed inside and outside of the planet’s orbit and 
in the corotation region are all included. Left and right panels re¬ 
fer, respectively, to the lower and higher gas densities (po = 10“^^ 
and 10“^^ gcm“^). The color scale indicates the logarithm of the 
number of particles. Inclinations are in degrees, radii are in cm 
and temperatures are in K. 

mals versus semi-major axis. The distributions include 
both icy and mixed-composition bodies, for po = 10“^^ 
(left) and 10“^^gcm“^ (right), initially deployed in all 
three disk regions mentioned above. As indicated by Fig¬ 
ure 7, orbital eccentricities can exceed 0.8 at both gas 
reference densities, but these values typically occur for 
^ 3 Up. Most of planetesimals produced via ablation 
of initially larger bodies have temperatures Tg > 150 K, 
which indicates that ablation is ongoing. In some in¬ 
stances, temperatures have settled at lower values, in¬ 
dicating that these planetesimals orbit in colder disk 
regions. At Tg > 200 K, bodies are consumed rela¬ 
tively rapidly. In fact, below the critical temperature 
{Tg < Ter), the timescale for ablation of a planetesimal 
is 


dMs 

^ 1 RsPs 

127rkBTg 

dt 

~ 3 Pv \ 

1 l^smu 


In the equation above, for mixed-composition planetes¬ 


imals, Py is the vapor pressure of ice and ps is the 
mean molecular weight of ice modified by the ice mass 
fraction. At Tg = 150 K, the vapor pressure of ice is 
^ 5 X 10“^dyne/cm^ and the timescale for ablation of 
a Rg = 100 km body would be in excess of 10^ years 
(assuming that Tg remains constant and neglecting re¬ 
condensation). But this timescale rapidly declines as 
Tg (and hence Py) rises, becoming 4-5 x 10^ years at 
200 K {Py « 1.6dyne/cm^) and ^ 10^ years at 210 K 
{Py ^ 7dyne/cm^). The distribution of Tg versus ag 
in Figure 10 would place the ice line at around 2.8 AU, 
where the gas temperature Tg = T^ is about 220 K (see 
Equation (3), although Tg needs not be in equilibrium 
with T^). However, global models of evolving disks (e.g., 
D’Alessio et al. 2005, see also D’Angelo & Marzari 2012) 
predict somewhat lower gas temperatures in those disk 
regions at times > 1 Myr. Therefore, the ice line can 
move inside 2 AU as long as the gas remains optically 
thick to stellar radiation. 

6.1. An Estimate of Collision Rates 

The rates of collisions among planetesimals can be de¬ 
rived from simple arguments. The average relative ve¬ 
locity between two bodies in a swarm is such that (e.g., 
Stewart & Kaula 1980) 

All) = . (34) 

Based on the two-body approximation, the cross section 
for collisions of a target planetesimal of radius Rg with 
those of radius Rj is 

5,.=7r(i?,+i?,y (^1 + ^), (35) 

where = 2GMg/Rg. The rate of collisions on the 
target body is then 

^ (36) 

3 

where Mj is the number density of solids, which involves 
the planetesimals’ surface density, and the swarm 
thickness, (agsinis). The averages {v‘^^f) and (agSinA), 
as a function of Og , are directly evaluated from the calcu¬ 
lations. In a swarm with equal number densities of 0.1- 
100 km bodies, the largest planetesimals have the high¬ 
est collision rates. In the circumstellar disk, at some 
distance from the planet’s orbit and for ^ lgcm“^, 
during the course of the calculations the number of col¬ 
lisions on any planetesimal would be negligible. The 
same arguments applied to planetocentric orbits in the 
circumplanetary disk (see Section 7) would yield a colli¬ 
sion rate of dNjdt ^ 10“^ GMp /(Es/lgcm“^) for 
Og ^ 0.15 Th- 

7. EVOLUTION OE PLANETESIMALS IN THE 
CIRCUMPLANETARY DISK 

Consider the reference frame {O'; x, z}, introduced 
in Section 5.1, whose origin is fixed to the planet. Based 
on relative positions and velocities, we derive the oscu¬ 
lating orbital elements in this frame of the planetesimals 
bound to the planet. Relative trajectories that are not 
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Figure 11. Planetocentric orbital elements, semi-major axis (left), eccentricity (center), and inclination (right), of icy and mixed- 
composition planetesimals bound to the planet. Top (bottom) panels refer to the density po = 10“^^ (10“^^) gcm“^. 


representable by osculating ellipses are disregarded, al¬ 
though they may belong to bodies that will accrete on the 
planet. Notice that we maintain the same notations as 
in previous sections, even though the orbits are relative 
to the planet. Since planetesimals are initially placed 
on orbits about the star well outside the Roche lobe 
region, these are all objects captured by the planet’s 
gravity. Capture is aided by dissipation of kinetic en¬ 
ergy through gas drag. Figure 11 shows semi-major axes 
(left), eccentricities (center), and inclinations (right) of 
these orbits, including both icy and mixed-composition 
bodies. The top/bottom panels refer to the lower/higher 
background density, po- The histograms include only 
those objects whose (planetocentric) orbit has a semi¬ 
major axis < 0.6 Rh- Inclination distributions indicate 
the presence of retrograde orbits {is > 90°, see also Fig¬ 
ure 12). These objects are among those initially moving 
in the corotation region and in the region exterior of the 
planet’s orbit, although the absence of retrograde objects 
originating from the inner disk region may be a result of 
small-number statistics. As explained below, since the 
orbital elements are based on osculating ellipses, not all 
planetesimals represented in Figure 11 are permanently 
captured (or accreted) by the planet. 

The histograms in the top row of Figure 11 include 
about 1.8 times as many bodies as the histograms in the 
bottom row. Yet, the distributions for the reference den¬ 
sity po = 10“^^ gcm“^ contain about 25% as much mass 
(see distributions of Rs in Figure 12). Most captured 
bodies have radii < 10 km, but most of the mass is car¬ 
ried by Rs ^ 100 km planetesimals. 

Average mass fractions on the order of several times 
10“^ are captured from the planetesimals initially placed 
in the inner disk and somewhat larger fractions, ^ 10“^, 
are captured from bodies initially deployed in the region 
beyond the planet’s orbit. Much smaller mass fractions 
are instead captured from the corotation region. The 
fractions are computed as the ratio of the mass of plan¬ 
etesimals moving inside r < 0.6 Rr to the total available 


mass in solids. Masses are averaged over the last ~ 50 
orbits of the planet. The mass in solids accreted by the 
planet is not taken into account. This fractional mass 
may be considered as an “equilibrium” mass between the 
supply of solids from the circumstellar disk and the mass 
loss due to ejection, ablation, break-up, and accretion 
of planetesimals in the circumplanetary disk. If results 
were rescaled so that the surface density of solids between 
0.77 ap and 0.82 and between 1.2 and 1.25 was 
1 g cm“^ at the end of the calculations, the average mass 
within 0.6 Rh of the planet would be r\j 10-3 Me. 

Planetesimals that orbit within ^ 20° of the equatorial 
plane are subject to an increasing drag force, as they ap¬ 
proach the planet, due to the augmenting gas density (see 
Figure 5). Equatorial gas rotation inside 0.1 Rr of the 
planet deviates only a few percent from Keplerian rota¬ 
tion Y^GMp/r, but the relative difference increases with 
distance, becoming 10% at r ^ 0.2 Rr and « 40% at 
f ~ 0.5 Rr. The radial velocity of the gas at the equator 
is much smaller in magnitude than the azimuthal veloc¬ 
ity. For near-circular orbits, gas density can be approxi¬ 
mated to a constant and Equation (32) may be applied, 
although the term in the square brackets of Equation (27) 
should depend on r in these cases. According to Equa¬ 
tion (32), the decay time of orbital semi-major axes due 
to aerodynamics drag at f ^ 0.5 Rr is on the order of 
a few times {Rs/as){ps/Pg), in units of the local orbital 
period around the planet. We recall that the length 
here represents the semi-major axis of the planetocentric 
orbit. The impact of gas drag on higher inclination orbits 
(40° ^ is ^ 140°) is probably somewhat smaller, since 
high densities are mostly encountered when these orbits 
are close the disk’s equatorial plane (see top panels of 
Figure 4). 

Figure 12 shows two-dimensional distributions analo¬ 
gous to those in Figure 10, but for planetesimals bound 
to the planet. Note that counts of 1, shown in Figure 12, 
are not visible in histograms of Figure 11, as they lie 
on the horizontal axis. The figure indicates that plan- 
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Figure 12. Distribution of icy and mixed-composition planetesi- 
mals bound to the planet versus semi-major axis and versus orbital 
eccentricity, orbital inclination, radius, and temperature, as indi¬ 
cated on the vertical axes. Osculating orbital elements are com¬ 
puted from positions and velocities relative to the planet. Left and 
right panels refer, respectively, to po = 10“^^ and 10“^^gcm“^. 
The color scale indicates the logarithm of the number of particles. 
Units are as in Figure 10. 

etesimal temperatures are typically Tg ^ 200-220 K (see 
also Figure 13). As argued above, the ablatiou timescale 
rapidly decreases at higher temperatures. A fiuer sam- 
pliug of the semi-major axis histograms reveals that uum- 
ber deusities start to decliue iuward of f ~ 0.05 i?H and 
~ 0.08 i?H foi* Po = 10“^^ aud 10“^^ gcm“^, respectively. 
Figure 13 shows that particle temperatures (blue circles) 
become cousisteutly lower thau gas temperature (solid 
hue) wheu substautial ablatiou begius, as expected from 
Equatiou (15) for slowly varyiug frictioual heatiug aud 
radiative gaiu/loss euergy terms. The gray circles be- 
loug to isothermal circumplauetary disk calculatious dis¬ 
cussed iu Sectiou 7.1. The large ablatiou rates cau eurich 
these circumplauetary disk regious with ice aud rock, lo¬ 
cally iucreasiug the solid-to-gas mass ratio. 

It is safe to assume that, over the course of the cal¬ 
culatious, solid material is practically ouly ablated ei¬ 
ther iuward of 2.8 AU or iu close proximity of the plauet. 
By separatiug these two coutributious, we estimate that 
^ 10“^MEyr“^ worth of ice aud silicates would be re¬ 
leased iu the gas close to the plauet, if the iuitial surface 
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Figure 13. Temperature versus semi-major axis of some planetes- 
imals bound to the planet (blue circles) for the reference gas density 
po = 10“^^ (top) and 10“^^gcm“^ (bottom). The light gray cir¬ 
cles indicate the temperature of planetesimals in the isothermal 
circumplanetary disk calculations (see Section 7.1). The solid lines 
represent the gas temperature Tg in Equation (31), also plotted in 
Figure 6. The dotted line is the gas temperature Tn in Equation (3) 
that, in the “isothermal” calculations, is applied everywhere in the 
disk. 

deusity of plauetesimals betweeu 0.77 aud 0.82 aud 
betweeu 1.2 aud 1.25 was 1 gcm“^. At this produc- 
tiou rate, the average metallicity of the circumplauetary 
disk would become ^ 0. 01(10 ^^gcm ^/po) in ^ 100 
plauet’s periods. There are ouly margiual differeuces iu 
the amouuts of material ablated from icy aud mixed- 
compositiou bodies (^ 10%, roughly cousisteut with the 
predictious from Equatiou (33) if equal temperature bod¬ 
ies are assumed) aud similar differeuces are obtaiued for 
the two values of pQ. 

Iu the calculatious, scatteriug of plauetesimals out 
of the Roche lobe occurs through iuteractiou with the 
plauet, followiug oue or more close eucouuters. Iu Eig- 
ures 11 aud 12, both low aud high ecceutricity orbits 
may eveutually lead to scatteriug eveuts. Of the objects 
plotted iu the upper pauels of Eigure 11, by followiug 
the subsequeut evolutiou of a sample, we estimate that 
the ratio of accreted to scattered (out of the plauet’s 
Roche lobe) objects is about 1/7. The Rg ^ 100 km 
plauetesimals, which are the least affected by gas drag, 
are those most easily ejected. Bodies are scattered both 
iuward aud outward of the plauet’s orbit. Iu this sam¬ 
ple, whose iuitial positious aud velocities are the same 
as those used to draw the distributious iu Eigure 11, the 
majority of scattered bodies with a startiug semi-major 
axis > 0.3 Rh have ecceutricities ^0.6. To give au idea 
of the spatial distributiou of scattered objects, Eigure 14 
shows the positious aloug the trajectories of 18 such bod¬ 
ies, of radius ^ 10 aud ^ 100 km (see figure captiou for 
details), extracted from said sample. Iu particular, the 
plauetesimals iu the figure are origiually placed exterior 
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Figure 14. Positions along 18 trajectories, in a reference frame 
fixed to the planet (Section 5.1), of icy planetesimals scattered 
out of the planet’s Roche lobe. The star is located at (x, y) = 
(—a,p,0). The reference density in the disk is po = 10“^^gcm“^. 
The trajectories are integrated for about 57 orbital periods of the 
planet. The color scale of the dots indicates the planetesimals 
radius (which increases from lighter to darker dots). Except for one 
planetesimal {Rs ~ 1km, scattered inward), about equal numbers 
of bodies have radii 10 and ~ 100 km. 
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Figure 15. Distributions of the radii of planetesimals accreted by 
the planet during the course of the calculations. The histograms 
include bodies of both compositions, initially placed in all three 
disk regions. The reference gas density is po = 10“^^ (top) and 
10“ii gcm“^ (bottom). The radii are in cm. Histograms are nor¬ 
malized so that the maximum is 1. 
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of the planet’s orbit, before being diverted toward the 
circumplanetary disk (but all trajectories illustrated in 
the figure start from within 0.6 Ru of the planet). 

Table 6 lists the accretion rates of planetesimals on the 
planet, estimated over the last 50-70 planet’s revolutions. 
The accretion rates are rescaled so that the average sur¬ 
face density of solids in the regions 0.77 < r/ap < 0.82, 
0.965 < r/ap < 1.035, and 1.2 < r/ap < 1.25 is 
lgcm“^. Contributions from each of the three regions 
are listed separately. There are relatively small differ¬ 
ences (< 10%) between the values obtained from cal¬ 
culations using different reference gas densities, po, and 
different material compositions. The accretion of solids 
arises almost entirely from the regions interior and exte¬ 
rior of the planet’s orbit, the former region contributing 
around 5% of the total. The corotation region provides 
only a minimal fraction of the solids’ accretion. Adding 
up the regional contributions and averaging out the re¬ 
sults, we have 

{Mj,)s = 2.8 X 10-5 ( ^ / ) Me yr'^, (37) 

\lgcm“^y 

where is the surface density of solids in the three 
radial regions mentioned above. This result applies for 
initial planetesimal populations with equal numbers of 
bodies per size bin. Additionally, as stressed in the pre¬ 
vious section, planetesimal-planetesimal interactions are 
neglected and, therefore. Equation (37) is valid as long 
as the effects of collisions and encounters among solids 
can be neglected, i.e., for low enough values of But 
in absence of a mechanism (like collisions and gravita¬ 


tional stirring) to replenish the exterior disk region with 
solids (resupply via gas drag would not be effective at 
the gas density levels considered here due to long orbital 
decay times, see Section 6), a total mass in solids of or¬ 
der 0.2 Eg/(I g cm“^) Me can be delivered to the planet. 
A somewhat smaller mass would be accessible from the 
interior disk, but over a much longer timescale. This 
amount of solids may represent only a relatively small 
addition to the heavy element content of the planet. 

The histograms in Figure 15 show the radii in cm 
of accreted planetesimals for po = 10“^^ (top) and 
10“^^gcm“^ (bottom). Histograms are rescaled so that 
the maximum count is 1. For the lower reference density 
case, there is an almost equal probability, within < 15%, 
of accreting bodies in the size range 1 km ^ Rs ^ 100 km. 
The difference increases for the higher reference density 
case, in which the probability of accreting Rg ^ 1 km 
bodies is about half that of accreting Rs ~ 100 km bod¬ 
ies. (As discussed below, one reason for the different 
accretion probability lies in the fact that smaller bodies 
are more likely to break up.) Consequently, our choice 
of starting the planetesimal populations with equal num¬ 
bers of objects in different size bins determines the result 
that the accretion rate in solids is basically supplied by 
the largest planetesimals. In fact, the accretion rate in 
Equation (37) is basically that of the largest planetesi¬ 
mals and Eg is the surface density of the largest bodies. 
Therefore, Figure 15 implies that, within factors of or¬ 
der unity. Equation (37) provides the accretion rates of a 
mono-size population of planetesimals (whose radius be¬ 
longs to the approximate range 1-100 km) with Eg being 
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Table 6 

Accretion Rates of Planetesimals on the Planet^ 


Icy Bodies Mixed Bodies 

Zone po = 10-^2 iq-ii ^ iq -12 ^q-ii 

Interior 1.2 x 10“® I.O x 10“® 1.2 x 10“® 1.5 x 10“® 

Corotat. 8.7 x 10“'^ 8.2 x 10“'^ 8.8 x 10“'^ 7.1 x 10“'^ 

Exterior 2.7 X 10“^ 2.6 x 10“^ 2.4 x 10"^ 2.7 x 10"^ 

^ In units of Mgyr”^ and scaled to Eg = Igcm“^. 


the surface density of that population. Figure 15 includes 
all planetesimals accreted during the course of the cal¬ 
culations, hence it contains some bias due to the choice 
of the initial distributions of planetesimals. For example, 
the probability of delivering to the circumplanetary disk, 
and thus of accreting, Rg ^ 0 . 1 km bodies declines over 
time because of the widening gap in solids of this size 
(see Figure 9). 

The initial semi-major axis, eccentricity, and inclina¬ 
tions of accreted planetesimals are plotted in Figure 16 
(see figure caption for details). The results indicate that 
planetesimals deployed in proximity of the edge closer 
to the planet’s orbit are more likely to be accreted than 
are more distant bodies (see left columns). For bod¬ 
ies released in the cororation region (not shown in the 
figure), the trend is opposite, as expected, due to the 
stability of tadpole orbits. The probability of accreting 
planetesimals deployed exterior of the planet’s orbit is 
fairly independent of the initial eccentricity, whereas (in 
the range 0-0.05) larger initial eccentricities favor accre¬ 
tion of bodies deployed interior of planet’s orbit (see cen¬ 
ter columns). Initially co-planar orbits lead to accretion 
more easily than do inclined orbits (see right columns). 
The smaller accretion rates provided by the interior disk 
(see Table 6 ) can be explained by observing that the peak 
number densities around the inner edge of the solids’ gap 
are at ^ 0.8 (see Figure 9), where the probability of 
accretion is relatively low (see Figure 16, left-even pan¬ 
els). 

As mentioned in the previous section, break-up of plan¬ 
etesimals may occur when they encounter the dense gas 
of the circumplanetary disk, inside 0.2 Ru of the planet 
(see Figure 5). Assuming that planetesimal fragments 
can quickly drift toward the planet, the mass of disrupted 
bodies would contribute to the accretion rate of solids 
on the planet. However, this addition would amount 
to ^ 0 . 01 % or ^ 0 . 1 % (depending on po? see above) 
of the values reported in Table 6 . This small contribu¬ 
tion largely depends on the fact that accretion of solids 
is dominated by Rs ~ 100 km planetesimals, which do 
not tend to break up. However, smaller bodies break up 
more easily. In fact, for po = 10“^^gcm“^, the mass of 
shattered planetesimals with radii 1 km ^ ^ 10 km 

is comparable to (and sometimes larger than) the ac¬ 
creted mass contributed by bodies of these sizes. Thus, 
if fragments ought to be considered as accreted mate¬ 
rial, the accretion rate in this size range (e.g., given by 
Equation (37) applied to a mono-size population) may be 
higher by a factor of up to a few. For the reference den¬ 
sity Po = 10 “^^gcm“^, break-up of planetesimals is less 
relevant, and it would contribute ^ 10 % to the accretion 
of Rs ^ 10 km bodies and even less to the accretion of 


smaller planetesimals. 

7.1. Isothermal Circumplanetary Disk Calculations 

As anticipated in previous sections, we also consider 
models that apply the temperature in Equation (3) 
also in the circumplanetary disk, which then becomes 
nearly isothermal with a temperature of ~ 120 K (see 
Eigure 6 at r ~ Rn and the dotted line in Eigure 13). 
A smaller number of bodies is released in these calcula¬ 
tions. Since the gas density is the same as in the models 
discussed above, differences may be expected especially 
in the thermal evolution and ablation of planetesimals 
moving in close proximity of the planet. Nonetheless, we 
find that there are not large differences between the two 
approaches. 

As above, average fractions of order 10 “^ are captured 
within 0.6 i?H of the planet from the available mass of 
planetesimals. More bodies are captured at the lower 
reference density, po, but a somewhat larger mass is re¬ 
tained at the higher value of po- Most of the captured 
bodies have radii < 10 km, but most mass is carried by 
the lOOkm-radius planetesimals. 

The gray circles in Eigure 13 show the planetesimal 
temperatures as a function of the planetocentric semi¬ 
major axis, while the dotted line represents the local 
gas temperature (see Equation (3)). In the calculations 
with reference density po = 10 “^^gcm“^, is gener¬ 
ally comparable with the planetesimals’ temperature ob¬ 
tained from the calculations discussed above (top panel, 
darker circles), which use the gas temperature in Equa¬ 
tion (31). Eor the higher reference density, the discrep¬ 
ancy is larger, but only when 0.03 Ru < Og ^ 0.2 Ru (see 
bottom panel). Body temperatures become comparable 
again (T^ > 200 K) when planetesimals get close to the 
planet (a^ < 0.03i^n), where ablation is most vigorous. 
It is possibly for this reason that the amounts of ablated 
material are similar in the two sets of calculations. 

The accretion rates of planetesimals on the planet are 
comparable to those given by Equation (37). Relative 
differences between isothermal and non-isothermal calcu¬ 
lations of accretion rates versus Rs are also small, < 20 %. 

The thermal distribution of the gas does not directly 
affect the fragmentation of planetesimals, as neither the 
dynamical pressure, Pdy, nor the material compressive 
strength, cTs^/I km/P^, is explicitly dependent on Tg 
{Tn) or Tg. Indirect effects may nonetheless occur, e.g., 
because of different ablation histories. Consistently with 
the results presented above, break-up of planetesimals 
occurs in the proximity of the planet, and mostly in the 
size range 1 km ^ Rs ^ 10 km. Overall, solid mate¬ 
rial made available through break-up would only con¬ 
tribute negligibly to the accreted mass due to the fact 
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Figure 16. Initial distributions of semi-major axis (left), eccentricity (center), and inclination (right) of accreted planetesimals for 
po = 10“^^ (upper pair of rows) and 10“^^ gcm“^ (lower pair of rows). The histograms refer to the initial distributions of both icy and 
mixed-composition bodies. Histograms are normalized so that the maximum is 1. 


that most mass is delivered to the planet in the form of 
Rs ^ 100 km bodies, which rarely break up. But this ef¬ 
fect depends on body size and gas density. As above, at 
po = 10“^^gcm“^ (but not at the lower po)? disruption 
of 1km ^ Rs ^ 10 km planetesimals would significantly 
contribute to the accretion of solids, if fragments effi¬ 
ciently accreted on the planet. 

8. DISCUSSION AND CONCLUSIONS 

The results presented here show how planetesimals of 
various sizes, initially orbiting in three narrow radial re¬ 
gions (Ar ^ 0.26-0.36 AU) around a star and in prox¬ 
imity of a Jupiter-mass planet, are scattered through 
the circumstellar disk and toward the planet. Scatter¬ 
ing is dominated by three-body interactions (star-planet- 
planetesimal) while gas drag typically operates as a per¬ 
turbing force. In some case, for Rs < 0.1km, gas drag 


does determine planetesimal dynamics. As expected, the 
evolution of mixed-composition (ice/quartz with a 60% 
ice mass fraction) and icy planetesimals is similar in most 
thermodynamical aspects. 

The surface temperature of planetesimals, (defined 
in Section 3.2), evolves toward an equilibrium value 
{dTs/dt ~ 0), which in absence of significant mass loss is 
such that (see Equation (15)) 

(3g) 

32crsB eg 

By using the approximation |vs — v^| ^^asUK/2 
with ^ H/r (which holds far from the planet’s 
orbit, see Section 5 and Appendix B.l), the second 
term on the right-hand side of Equation (38) becomes 
g/2)®C'D/(4crsB)(/Cg/es)(GM*/as)^/2. In the calcula- 
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tions described here, this term is usually small compared 
to Tg and thus ^ Tg. This conclusion is supported 
by the calculations. Vigorous mass loss tends to lower 
Ts relative to Tg (see Figure 13). The inverse of the 
timescale required to reach equililDrium, neglecting both 
frictional heating and cooling via latent heat release, is 


1 

dTs 

(JsbTs 

1-f^V 

Ts 

dt 

^ Ss PsCs 



For a given deviation from equilibrium (T^), this 
timescale is oc Sg/T^ oc Xg/T^ (see Equation (14)) and 
increases with approaching equilibrium. Some numerical 
examples on the evolution of Tg toward equilibrium are 
shown in Appendix B. 2 . The situation is more complex 
for eccentric orbit bodies, as they experience a varying 
gas temperature along their orbit. Nonetheless, an equi¬ 
librium temperature can be reached if \dTg/dt\/Tg ^ Ok- 

Results concerning the distributions of solids can be 
rescaled to an arbitrary surface density (in the radial 
regions of initial deployment), Eg, provided that colli¬ 
sions and encounters among bodies do not significantly 
alter their dynamics, i.e., that remains relatively 
low, as it may be the case at the late epochs of giant 
planet formation (see Section 6.1). The ejection rate 
of mass out of the disk domain ( 2 AU ^ r < 20 AU) 
is ^ 2 X 10“^ (Eg/l gcm“^) Me yr“^. In reality, most 
of these bodies are on orbits bound to the star when 
they cross the boundaries, and “ejection” generally clas¬ 
sifies orbits whose perihelia (aphelia) are inside (outside) 
of rmn (^mx)- Planctcsimals are scattered from the in¬ 
terior to the exterior of the planet’s orbit at a rate of 
^ 5 X 10“^ (E 5 /I gcm“^) Me yr“^, and in the opposite 
direction at a rate of ^ 2 x 10“^ (E^/l g cm“^) Me yr“^. 
These rates refer to scattered objects that have ellip¬ 
tical orbits about the star. Arguably, these scattering 
rates may only apply for a limited period of time, if 
solids are not replenished via collisions or some other 
mechanism (e.g., gravitational stirring or scattering by 
other planets). Therefore, the minimum masses that 
can be scattered out of boundaries, inside and out¬ 
side the planet’s orbit are r\j 0.25 (Ss/l gem ^)Me, 
~ 0.3(Ss/lgcm“2)ME, and ~ 0.1 (Ss/lgcm”^) Me, 
respectively (neglecting contributions from the mass in 
the corotation region). 

For the disk temperatures applied here, both icy and 
mixed-composition bodies would be ablated inside r ^ 
2.8 AU {Tg = Tji ^ 220 K). However, disk models (e..g, 
D’Alessio et al. 2005) suggest lower temperatures after 
a few to several million years, hence ice-rich planetesi- 
mals may survive at radii r < 2AU. Regardless, scat¬ 
tering by a Jupiter-mass planet provides an important 
source of hydrated planetesimals to inner disk regions. 
A mass equal to the current mass of the main aster¬ 
oid belt (Krasinsky et al. 2002) would be delivered in 
^ 50 (1 gcm“^/Es) yr. These bodies would still orbit 
in a relatively dense gas, but the orbital decay time 
of Rg > 1 km planetesimals around 2.5 AU would be 
> 5 X 10^ local orbital periods for po ~ 10“^^gcm“^. 
Orbital eccentricities and inclinations would be damped 
on timescales shorter by factors of ^ 200 . 

The planetesimals orbiting in the corotation region 
ap ± Rh are removed at an average rate of ^ 6 x 


10“^ (E 5 /I gcm“^) Me yr“^. However, the tadpole or¬ 
bits around the L 4 and L 5 points are very stable, due to 
low gas densities (see Figure 2 ), and number densities do 
not drop around these points. Since the rate of capture 
in the corotation region appears to be much smaller than 
the removal rate, the local density may bear information 
about dynamical and physical conditions at earlier times, 
before the giant planet acquires its massive envelope (see 
also Peale 1993). 

The steep radial pressure gradient induced by the 
planet at the edges of the gap in the gas density profile 
(as function of r, see Figure 3) can partially prevent small 
planetesimals {Rg ^ 0 . 1 km) from crossing the planet’s 
orbit, but this effect reduces as gas dissipates (see Fig¬ 
ure 9). The size range most affected is determined by 
the strength of the drag acceleration (Equation (9)), and 
hence by dPg/dr at the gap’s outer edge. A condition 
for gap formation derived from the balance of viscous and 
tidal torques (D’Angelo & Marzari 2012 , and references 
therein) is 



where A = max (M, Rr)- The value of the left-hand side 
of this inequality is ~ 5.6 for the parameters adopted 
here. A similar number is obtained for a Saturn-mass 
planet and somewhat smaller values for H/r and 
which suggests that a reduction in the inward flux of 
Rg ^ 0.1 km (and smaller) planetesimals may begin prior 
to reaching the current mass of Jupiter. 

Experiments conducted on 1 cm < < 10 m bod¬ 

ies, initially released exterior of the planet’s orbit ( 1.2 < 
agjap < 1.25), show that these particles remain segre¬ 
gated. After ^ 600 orbital periods of the planet, results 
indicate that the amount of solids delivered to the inte¬ 
rior disk is negligible (for po = 10 “^^ gcm“^) or virtually 
zero (for po = 10“^^gcm“^). The radial position of the 
swarm’s inner edge is between r ^ 1.27 and 1.4 
for 1cm ^ Rs ^ 10 cm, and between r ^ 1.36 and 
1.4 Up for 1 m < Rg < 10 m (the position also depends on 
po). Eor po = 10“^^gcm“^, particles of 1cm in radius 
are halted at r ~ 1.18 a^, close to the peak of super- 
Keplerian rotation in Eigure 3. Segregation also leads to 
negligible fluxes of these bodies toward the circumplan- 
etary disk. When po = 10“^^gcm“^, we do find that 
10 m-bodies can be scattered toward the planet and the 
inner disk when ag ~ 1.2 and eg ^ ig ^ 0. However, 
this scattering event lasts only briefly at the beginning 
of the calculation, before the swarm recedes. Therefore, 
this is likely a transient effect induced by the choice of 
the initial distributions. Nonetheless, if ^ 10 m bod¬ 
ies are produced via collisional comminution of planetes¬ 
imals around r = 1.2 a^, part of them may be scattered 
inward. 

Planetesimals can be captured in the circumplanetary 
disk with a wide range of (planetocentric) orbits, includ¬ 
ing retrograde ones (see Eigures 11 and 12 ) as also found 
by other recent studies (e.g., Eujita et al. 2013; Tanigawa 
et al. 2014). The ensemble of bodies with retrograde or¬ 
bits comprises planetesimals with radii 0 . 1 km ^ Rs ^ 
100 km, although most of those coming from the exterior 
disk have Rg > 10 km. Capture of planetesimals provides 
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the circumplanetary disk with a time-averaged solids’ 
reservoir of ^ 10 “^ (E^/l g cm“^) Me, which may be 
considered as a balance between the external supply and 
the loss due to ejection, ablation, break-up, and accretion 
on the planet. This amount of solids would account for 
a relatively low surface density (although planetesimals 
are continuously supplied), ^ 0.3 which may indicate 
relatively long times, ^ 10 ^ (1 gcm^^/Eg) local orbital 
periods about the planet, for the formation of ^ 10 ^ km- 
radius satellites. Applied to the Galilean satellites, these 
formation times appear compatible with Callisto, which 
is partially undifferentiated (Stevenson et al. 1986), and 
may suggest post-formation differentiation of the inner 
three satellites (Schubert et al. 2004). Type I migra¬ 
tion due to tidal interactions (e.g., Tanaka et al. 2002) 
with the thick circumplanetary disks considered here (see 
Figures 4 and 5) would lead to timescales for the orbital 
decay, at f ^ 0.04 i?H, of ^ 10 ^ (10“^^ g cm“^/po) local 
orbital periods, longer than formation timescales. 

Sustained ablation close to the planet (f < 0.1 i^n) 
releases heavy elements in the gas at a rate of ^ 
10 “^ (E 5 /Igcm“^) Me yr“^, which is large enough to 
significantly alter the gas metallicity over relatively short 
timescales, possibly leading to a dust laden system. 
Disruption of planetesimals may also contribute to the 
solids’ content of the circumplanetary disk. 

Equation (37) approximates the accretion of solids on 
the planet supplied by a mono-size swarm of planetesi¬ 
mals, with radius in the range from ^ 1 km to ^ 100 km, 
where is the solids’ surface interior and exterior of the 
planet’s orbit. Figure 16 indicates that the efficiency of 
accretion of planetesimals declines with increasing dis¬ 
tance from the planet’s orbit. If the edges of the gap 
in the solids’ distribution are eroded, because of lack of 
supply or because they recede due to gas drag torques, 
the accretion rate is expected to decrease. Probably, late 
accretion of solids only represents a minor addition to the 
heavy element content of a giant planet (unless E^ is still 
quite large). 

We estimate the mean accretion energy per unit mass, 
(AEacc/^Afg), delivered to the planet by accreted plan¬ 
etesimals during the course of the calculations. The 
quantity AEacc contains both kinetic and gravitational 
energy. We assume that all this energy is delivered close 
to the planet surface. The energy per unit time produced 
by accretion of solids is then 

{AE^JAM,){Mp), ~ 10-5 ( Lq. (41) 

\ 1 g cm ^ y 

This accretion power can be compared to the planet’s 
luminosity due to envelope contraction, between ^ 
10“^ I/Q and ^ 10“^ Lq (Lissauer et al. 2009). In reality. 
Equation (41) gives only a lower limit to the accretion 
power since additional energy is released as planetesi¬ 
mals sink into the planet. In fact, while ice dissolves 
at relatively shallow depths, where the temperature is 
< Ter ^ 650 K, rock (which makes 40% of the mass 
of mixed-composition planetesimals) can sink to much 
deeper layers on account of the higher critical tempera¬ 
ture (Ter = 4500K for quartz). 

In some instances (1km ^ Rs ^ 10 km and po ^ 
10 “^^ gcm“^), planetesimal break-up in the circumplan¬ 
etary disk produces significant amounts of solids that 


may increase {Mp)s by factors of up to a few, if debris 
is accreted before being completely ablated (the dissolu¬ 
tion timescale is oc see Equation (33)). However, if 
break-up produces large (^ 0 . 1 km) fragments, ablation 
appears to dominate over accretion by a large margin: 
while ^ 50% of their mass is ablated, only ^ 0.01% is 
accreted. To examine more in detail the fate of smaller 
fragments, we present tests that use 1 cm < < 10 m 

bodies as a proxy. These are released on circular orbits 
around the planet, between r ^ 0.1 Tr and ^ 0.6 Tr, 
at the disk’s equator. We find zero or negligible accre¬ 
tion. For radii 1 m < Tg < 10 m, planetesimals migrate 
inward very quickly, but not enough to overcome abla¬ 
tion. Essentially, they are all ablated. At the reference 
density po ^ 10 “^^gcm“^, disruption of planetesimals 
releases less mass, but fragments can still be produced 
via collisional comminution. The same tests reveal sim¬ 
ilar conclusions. Bodies with radii Im < Rg < 10 m are 
almost entirely ablated without any significant amount 
being accreted. For comparison, Rs ^ 0.1km planetes¬ 
imals shed in the gas via ablation ^ 30 times the mass 
they deliver to the planet via accretion. 

Solids in the range 1 cm ^ < 10 cm, for both val¬ 

ues of po, are also much more prone to ablation than 
they are to accretion, if they move toward the planet. 
In the tests, none of these particles is accreted. In ei¬ 
ther case, the mass that is not ablated remains beyond 
f 0.1 Tr for the duration of the calculations, possibly 
because these small solids are more efficiently coupled to 
the gas than are larger particles and the gas radial veloc¬ 
ity T'Vgjr (where Vg is relative to the planet) is positive 
at the equator, i.e., directed away from the planet (see 
also Tanigawa et al. 2012). The conclusion is that small 
fragments resulting from disruption should not signifi¬ 
cantly contribute to accretion, but should rather con¬ 
tribute to the local reservoir of solids and to enriching 
the gas with heavy elements. 

The main limitation of this study is the lack of 
planetesimal-planetesimal interactions, especially in the 
circumplanetary disk (see Section 6.1), which could affect 
the distribution of solids but which allows us to rescale 
the outcomes of the calculations to different values of 
the initial surface density of solids. The relatively short 
time span covered by the models is also a limiting factor. 
Another limitation is obviously the “discrete” approach, 
i.e., that of treating each particle as an individual body, 
which prevents from dealing with more realistic swarms 
of planetesimals, in terms of both number densities and 
size distributions. However, this approach allows us to 
model the evolution of the thermodynamical properties 
of single planetesimals at levels of detail not accessible to 
other, e.g., statistical or hybrid, approaches. Therefore, 
the method applied here can complement other tech¬ 
niques by providing detailed information on restricted 
populations of planetesimals at selected epochs of evolu¬ 
tion. 
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APPENDIX 

A. THE DRAG COEEEICIENT 

Melosh & Goldin (2008) performed an extensive study 
of the literature on existing gas drag experiments (see, 
e.g., Walsh 1976, and references therein). They derived 
an expression for the drag coefficient in Equation ( 8 ) 
and (9), as a function of the Mach number A4 (see Equa¬ 
tion (10)), the Reynolds number 7Z (see Equation (11)), 
and their ratio 1C = MjlZ. The function is continuous 
and extends over the entire (plausible) ranges of M. and 
7^. 

Let us introduce the adiabatic gas sound speed 


c 


9 ~ 



llgiriYL ’ 


(Al) 


The ratio 1C of the Mach number to the Reynolds number 
can also be written as 


^ 32V^(4) 

It is important to note that 1C is proportional to the 
Knudsen number, which is defined as the ratio between 
the mean-free path of a gas atom/molecule and the par¬ 
ticle diameter. Therefore, 1C may be regarded as a mod¬ 
ified Knudsen number. The proportionality factor de¬ 
pends on the form adopted for the dynamical viscosity 
rjg (see PPR 88 ). In our case. Equation (A4) yields a pro¬ 
portionality factor equal to (16/5)^7p/(27r) 1.28.^/^. 

In the derivation of Melosh & Goldin (2008), the drag 
coefficient is written as 


= 2 + (Cs - 2) + Ce 

where 



(A7) 

(AS) 


where 7 ^ is the adiabatic index of the gas (the isothermal 
sound speed is obtained for 7 ^ = 1 ) and the mean ther¬ 
mal velocity of the gas (e.g., Mihalas & Weibel Mihalas 
1999) 



TT 


(A 2 ) 


which is the equivalent of Equation (18) for the gas con¬ 
stituents (atoms and/or molecules). Let ns now define 
the magnitude of the relative velocity between the gas 
and a solid particle as = |v^ — v^j, then the (relative) 
Mach number can be written as 

The definition of the relative Reynolds number in 
Equation (11) involves the dynamical molecular viscos¬ 
ity of the gas. If interactions among gas atoms/molecules 
can be described as collisions between two rigid elastic 
spheres, an approximation of the dynamical molecular 
viscosity is (Mihalas & Weibel Mihalas 1999) 


where mn is the hydrogen mass and dn is the typi¬ 
cal diameter of the gas constituents. This length is 
dn = 2.71 X 10“^ cm for hydrogen molecules and 2.15 x 
10 “^ cm for helium (Haynes 2011). Although the inter¬ 
action model based on the rigid sphere representation of 
gas constituents, which interact only upon “contact”, is 
rather simple (see discussion in Mihalas & Weibel Miha¬ 
las 1999), Equation (A4) agrees within 25% with molec¬ 
ular hydrogen viscosity data (and 20 % with helium) in 
the temperature range from 100 to 600K (Haynes 2011). 

By substituting Equations (A4) and (A3) into Equa¬ 
tion ( 11 ), one finds that the (relative) Reynolds number 
can be cast into the following form 


7^ = 


_ 32^ / d^ 
\mu 


1^9 


PgRs^Ci.. 


(A5) 


and the auxiliary function G{1Z) is such that 


logG 


2.5 (7^/312)^" 

TT(W31^‘ 


(A9) 


The constants pi in Equation (A7) and p 2 in Equa¬ 
tion (A9) are pi = 3.07 and p 2 = 0.6688. The func¬ 
tion G{Tl) takes limiting values of 1 , for 7^ ^ 0 , and of 
lO^-^ ^ 316.23, for 7^ ^ 00 . 

Eor /C ^ 1 , when the particle size is much smaller 
than the mean-free path of the gas constituents, a regime 
referred to as free-molecular flow, the drag coefficient 
takes the value 

Cd Ce + 2 , (AlO) 

which, for Mach numbers ^ 1, becomes (4.6 + 

\.7^Ts/Tg)/ as in the Epstein regime (e.g., 

Whipple 1973; Weidenschilling 1977). We stress here 
that Equation (A 8 ) ought to be regarded as an exten¬ 
sion of the Epstein drag coefficient (see Hood & Horanyi 
1991 and discussion in Liffman & Toscano 2000). In fact, 
the form of Epstein coefficient typically adopted in the 
literaure, (8/3)^ 8 / {'K^g)/M (e.g., Supulver & Lin 2000; 
Chiang & Youdin 2010, and references therein), only ap¬ 
plies when there is specular reflection of the gas con¬ 
stituents impinging on the particle (see Epstein 1924, for 
details), which corresponds to assuming = 0 in the 
limiting expression above. Epstein (1924, Part I, Sec¬ 
tion 7) also argued that, for /C 1, particles should be 
considered as perfect thermal conductors, i.e., = T^, 

and the drag coefficient is then ( 8 / 3 + 7 r/ 3 )y^ 87 ( 7 n^/Al, 
in agreement with the limiting expression above. At large 
Mach numbers, the right-hand side of Equation (AlO) 
has asymptotic behavior 2 + (1.1 /M)^JTs/{'^gTg) (e.g., 
Whipple 1950; Baker 1959). 

Eor /C ^ 1 , as happens in the continuum and incom¬ 
pressible (M 1 ) flow regimes, the first exponential in 
Equation (A7) tends to 1 while the second exponential 
tends to 0. Therefore, the drag coefficient takes the value 


Gd 


/C<1 


Gs. 
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(All) 
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Figure 17. Top-left: the drag coefficient, Cd, given in Equation (A7), versus the Reynolds number, IZ (bottom axis), and K = J^jlZ 
(top axis), for Tg = Tg and a Mach number Ai = 0.01. The dotted line is the limit for /C ^ 1, given in Equation (AlO). The dashed line 
is the limit for /C 1, given in Equation (A12). Top-right: Cd as a function of the Mach number for particle of different radius: 10^ cm 
(solid line), 10^ cm (long-dashed line), and 10^ cm (short-dashed line). Bottom-left: The solid line is Cd from Equation (A7) while the 
dashed line is the drag coefficient used by Whipple (1973) and Weidenschilling (1977) (Tg is set to zero in Equation (A7), as discussed 
in the text). The inset shows a comparison with the drag coefficient used, among others, by Stalder & Zurich (1951); Probstein (1968); 
Hood & Horanyi (1991); Tedeschi et al. (1999); Liffman & Toscano (2000) (solid circles), which applies for /C ^ 1. Bottom-right: as in the 
bottom-left panel, but the dashed line is the coefficient for IZ > 1 and Ai > 1 used by PPR88. 


In this study, for Cs, we use a formula suggested by for free-molecular flows of Stalder & Zurich (1951, their 
Brown & Lawler (2003) Equations (A15) and (A17)). 


Cs 


24 

Tl 


(1 -f 0.1577^°) 


0.40777 

77-1-8710’ 


(A12) 


B. TESTS ON SOLUTIONS OE THE PARTICLE 
EVOLUTION 


in which the constant in the power of 7^ is ps = 0.681. For 
7^ < 1, Equation (A12) becomes the classical Stokes drag 
law Cs ~ 24/7^ (e.g., Whipple 1973; Weidenschilling 
1977; Brown & Lawler 2003), whereas, for 7^ ^ 1, we 
have that Cs ~ 0.407, sometimes referred to as the New¬ 
tonian drag coefficient (e.g., Whipple 1973). For non- 
spherical shapes, e.g., a cube or a short cylinder, this 
asymptotic value would be more than twice as large. 

In the top panels of Figure 17, we plot the drag co¬ 
efficient in Equation (A7) versus the Reynolds num¬ 
ber (Equation (A5)) and the modified Knudsen number 
(Equation (A6)), and also display the two limiting cases 
in Equations (AlO) and (A12). In the right panel, Cd is 
plotted for three different particle radii versus the Mach 
number. In the bottom panels, we make comparisons 
with drag coefficients used in previous studies (see figure 
caption for details), including the widely used coefficient 


In this Appendix, we present tests of the ordinary dif¬ 
ferential equation solver applied to the system of eight 
Equations (26), (4), (5), (6), (15), and (19) or (23). In 
order to make comparisons with compact analytic solu¬ 
tions, in the various tests we solve a reduced system and 
discuss separately dynamical problems (Equations (26) 
(4), (5), and (6)) in Appendix B.l and thermodynami¬ 
cal problems (Equations (15) and (19) or (23)) in Ap¬ 
pendix B.2. 

B.l. Dynamics Tests 

The solver is first tested against standard two-body 
problems, in which the particle orbits the star. Orbital 
energy and angular momentum are expected to be con¬ 
served in these problems and the extent to which this 
requirement is fulfilled provides an indication of accu¬ 
racy. 
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Figure 18. Variations of orbital angular momentum, semi-major axis and eccentricity in two-body problems with initial eccentricites 
0 (left), 0.5 (center), and 0.99 (right). In the left panel, to separate the three curves, Aag and Ae^ are shifted by an amount equal to 
5 X 10“^^. In the center and right panels, ALg is zero within machine precision, hence it is shifted by 2 x 10“^^ to appear in the plot. 


Let US indicate with the stellar mass, with Mg the 
particle mass, and with and the particle’s semi¬ 
major axis and eccentricity. In a two-body problem, the 
orbital energy and angular momentum per unit mass are, 
respectively. Eg = —G {M^ + M^) /{2as) and 

L, = v/G(M* + M,)a,(l-e2), (Bl) 

where G is the gravitational constant (here Lg should 
not be confused with the specific vaporization energy). 
Conservation of energy and angular momentum trans¬ 
lates into constancy of and or of and Lg. In fact, 
taking the differential of Equation (Bl) and dividing by 
we have 

dL, ^Idas ( el \ de^ 
i, 2 a, e,’ ^ > 


which connects the relative variations of ALg/Lg, 
Aag/ag, and Aegjcg (e.g., Beutler 2005). 

Experiments indicate that an advantage of integrating 
the equation of motion in terms spherical polar coordi¬ 
nates and angular momenta, in place of the usual carte¬ 
sian positions and velocities, is a substantial improve¬ 
ment in conservation of angular momenta and, typically, 
of energy. In the calculations reported in Eigure 18, we 
consider orbits with eccentricities = 0 (left), 0.5 (cen¬ 
ter), and 0.99 (right). In all cases Lg is conserved to 
machine precision, whereas the expected error in energy 
for the most eccentric orbit is one part in 10^ over a pe¬ 
riod of 1 Gyr, as the asymptotic error is linear in time 
in that case (e.g., Calvo & Sanz-Serna 1993). If neces¬ 
sary, better conservtion can be obtained by constraining 
the internal time step of the solver at the expense of an 
increased run time. 

Another test we discuss is a circular restricted three- 
body problem, constituted by two massive bodies, whose 
masses are Mi and M 2 , and a massless particle. All bod¬ 
ies orbit in the same plane and the radius of the massive 
bodies’ orbit is a. The Jacobi’s integral of motion for 
such system is (e.g., Murray & Dermott 2000) 


Gj = 


G(Mi+M2) 


r2 + 2 


GMi 

ri 


^ GM2\ 

H-I - ' 

r2 ) 


(B3) 

where r and v are the distance and velocity of the particle 
relative to the center of mass of the massive bodies, and 
ri and r 2 are the distances relative to these bodies. 

In Eigure 19 (left), we set Mi + M 2 = 1 and M 2 /Mi = 


0.001. The Jacobi’s integral is plotted as a function of 
the orbital period of M 2 around Mi, for particles on 
tadpole (thicker curve), horseshoe, and circulating (thin¬ 
ner curve) orbits (see figure caption for details). While 
the error in the circulating orbit test displays a typical 
asymptotic linear behavior (usually due to truncation er¬ 
rors in the algorithm), no systematic errors appear in the 
solutions for the tadpole and horseshoe orbits. 

In order to test the solver in the presence of drag, we 
follow the approach of Peale (1993). Consider a particle 
orbiting a star in gaseous disk. The rate of change of the 
particle’s orbital energy is equal to the work done on it in 
the inertial frame, that is MgdEgjdt = FD'^/g^ where 
is the drag force given by Equation (8). Differentiating 
the specific orbital energy Eg (see above), we have 

G {M^ ^ Mg)~\ / dag\ 3 Gd f Pg \ \ 1 

[—sg—J 

X (vg-v, - |v,|^) . (B4) 

Eor the sake of simplicity, the drag coefficient is taken to 
be constant, the disk’s gas velocity is approximated as 
sub-Keplerian (due to support provided by the pressure 
gradient) with no radial component, and jv^j agftKi 
where = G(M^ + Ms)/a^. Hence, we have that 

I Vs — v^l = Us fix (1 — a/I — C^). The quantity ^ is con¬ 
nected to the gradients of temperature and surface den¬ 
sity of the disk’s gas, as well as to the disk’s local thick¬ 
ness, H/r (see, e.g., Peale 1993; Takeuchi & Lin 2002; 
Tanaka et al. 2002), and is assumed to be constant. Un¬ 
der typical disk conditions, one finds that ^ ^ H/r (see 
also Equation (27)). If we indicate with ao and Uq the 
initial values of ag and Uk, Equation (B4) can be written 
as 


d 

dt 





(B5) 


in which 1/r = {3/4:)CD{Pgo/Ps){ao/Iis)^o and pg = 

PgO {cLo/daf with b > 0. The solution of Equation (B5) 

is 


Qjg 

ao 


1 + 25 


(i - 


2/(l+26) 


(B6) 

Our working assumptions imply that \das/dt\ <C 
asflK(l — yi — which, by using Equation (B5), be- 
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Figure 19. Left. Variation of the Jacobi’s integral (Equation (B3)) versus the orbital period of the massive bodies. The three curves refer 
to particles on different types of orbits: tadpole (thicker line) horseshoe, and circulating (thinner line). Center. Difference |Aas| between 
calculated and predicted positions of a particle subject to gas drag versus time in units of 27r/Oo. The predicted position is given by 
Equation (B6) with parameter b = 0. The three curves refer to different values of the constant Clor, as indicated in the legend (r is defined 
after Equation (B5)). Right. Difference \Avs\ between calculated and predicted (Equation (BIO)) velocities of a free-falling particle. The 
curves correspond to three values of the transient time r (defined after Equation (B9)) in units of l/12o, as indicated in the legend. 


comes (ao/r)(l - - i‘^){ao/aaf or 

equivalently 


3 

4 


Cd 



(t) (‘-v^)(S)* 


If ^ 1 and 0 < 6 < 1, the inequality (B7) becomes 

{‘i/S)CD{Pg/Ps){ao/R,)e < 1 - 

In the center panel of Figure 19, we show results for 
the orbital evolution of a particle initially moving on a 
circular orbit and subject to gas drag. The difference 
Atts between the calculated position and that predicted 
by Equation (B6) is illustrated for three values of ftor, 
assuming a radially constant gas density pg (i.e., 6 = 0). 

We also present a test conducted on a classical free-fall 
problem. Consider a particle at some height above the 
equatorial plane of the disk and suppose that it is sub¬ 
ject to a constant gravitational acceleration, directed 
toward the equatorial plane, and to gas drag (Equa¬ 
tion (9)). Moreover, suppose that the particle’s velocity, 
Vg, is perpendicular to the disk’s equatorial plane and 
that there is no vertical motion of the gas. The scalar 
acceleration of the particle is then 


_ _ _ 3 /Cp pg\ 

dt ^ ^\Rsps) 


\Vs\Vs‘ 


(B8) 


If Us >0, the acceleration is always negative, and even¬ 
tually the velocity becomes first zero and then negative 
{dvs/dt < 0 if Us = 0). If Us < 0, Equation (B8) can be 
written as 


dvs 

dt 



(B9) 


where l/igr'^) = {3/2){CD/R3){pg/Pa) and r is a 

timescale. For Vg < —2^r, the acceleration is positive 
and negative otherwise. Thus, the particle will always 
approach the velocity — 2^r, which is referred to as ter¬ 
minal or asymptotic velocity of the free-fall problem. The 
solution to Equation (B9) is 


Vs 


2gr 


\1± ) ‘ 


(BIO) 


Quantity 5 is a positive integration constant determined 
through the initial condition. The top (bottom) signs in 
front of B apply if u^ is smaller (larger) than {2gr)‘^. In 


the limit t ^ oc, either solution tends to the terminal 
velocity. 

In the right panel of Figure 19, we plot the difference 
Aus between the calculated free-fall velocity of a particle 
(with zero initial velocity) and that predicted by Equa¬ 
tion (BIO) as a function of the normalized time t/r. A 
velocity within 1% of the terminal velocity is attained for 
t/r > 5. 


B.2. Thermodynamics Tests 

In this section, we test the numerical solution of re¬ 
duced forms of Equation (12) against analytical solu¬ 
tions. We shall assume that heating and cooling pro¬ 
cesses affect the entire volume of the particle, which thus 
has a uniform temperature throughout. This assumption 
basically implies that the thermal conductivity of the 
body. As, tends to infinity (see discussion in Section 3.2). 
We adopt this simplified approach here, instead of solving 
Equation (15), because it helps in searching for analytical 
solutions of the reduced equations. 

Consider the situation in which the particle mass, M^, 
is constant and the temperature of the gas, T^, is always 
equal to T^, the particle temperature. Hence, if the par¬ 
ticle moves through gas around a star, it is constantly 
heated via friction so that its temperature changes in 
time according to 


d^ 

dt 


/ CpPg \ 

32 \RsPaCsJ 


(Bll) 


All quantities in parenthesis on the right-hand side are 
taken as constants. Following one of the problems in Sec¬ 
tion B.2, we assume that the gas is partially supported 
by pressure and has no radial velocity component, then 
I Vs — I = asf2K(l — a/I — ^^), where the azimuthal 
velocity of the particle is equal to If the radial 

position of the particle is given by Equation (B6) with 
6 = 0, Equation (Bll) becomes 


dTs _ 3 / CpPg \ /aof2o\/CT V 
dt ~ 32 \R^p,Cs) V i-ct/2 ) 


(B12) 


where r is defined beneath Equation (B5) and c = (1 — 
r. As above, oq and are the initial values 
of Og and Note that condition (B7) applies since we 
are using Equation (B6). The solution of Equation (B12) 
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Figure 20. Left. Difference |ATs| between calculated and predicted temperature of planetesimals, divided by the predicted temperature. 
Planetesimals are only heated via gas friction and cannot cool. The predicted temperature is given by Equation (B13). The curves refer to 
different values of the constant as indicated, and r is defined after Equation (B5). See text for further details. Center. Temperature 

evolution of planetesimals heated by the gas radiation field and losing energy via radiative cooling. The gas temperature is fixed at 
Tg = 150 K. The initial temperature of the planetesimals is 100 K and 200 K and their radius in km is indicated in the legend. The lower 
panel shows the normalized difference between the numerical and the analytical solution (Equation (B15)). Right. Thermal evolution of 
planetesimals that lose energy due to ablation. The temperature is shown on the top while the normalized difference between numerical 
and analytical temperature (Equation (B19)) is shown on the bottom. In all panels, the time units are 27t/CIq. 


is 


whereas for the function is 


T. = r,(o) + | 


/ CpPg \ 
\ RsPsCs ) 




1 1 

c(l —ct/2)2 c 


(B13) 


A Taylor expansion around t = 0 of the function in 
square brackets on the right-hand side gives t. The tem¬ 
perature diverges as t 2/c, i.e., as ^ 0, since the 
relative velocity jv^ — Wg\ diverges. 

We solve Equation (B12) for different values of the con¬ 
stant Qqt and initial temperature of 100 K. In Figure 20 
(left), we plot the difference |ATs|, between numerical 
and analytic solution (Equation (B13)), divided by the 
analytic solution. The rise in temperature is limited to 
a few degrees in case with longest r and to over 50 K in 
the opposite case. 

Consider another situation in which the particle has 
again a constant mass but its temperature differs from 
the gas temperature. Hence, the particle experiences 
heating by gas-emitted photons and cooling via back- 
body emission. Suppose also that the frictional heating 
is negligible or, otherwise stated, that the quantity in 
parenthesis on the right-hand side of Equation (Bll) is 
vanishingly small. Thus, the temperature variation of 
the particle is governed by 


dt 



(B14) 


Equation (B14) implies that the particle temperature 
evolves towards Tg, which is assumed to be constant. All 
quantities in first set of parenthesis on the right-hand 
side are also supposed to be constants. A solution to the 
equation is 


I{Ts) 


1 

4T3 
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In 


\Tg-Tj 


+ 2 arctan 



• (B17) 


Equation (B15) defines implicitly Tg as a function of 
time. 

Equation (B14) is solved numerically for three differ¬ 
ent radii (see figure caption) of planetesimals, orbiting a 
star with a period 27r/^4o- Two values of the initial plan- 
etesimal temperatures are applied: 100 K and 200 K, so 
that the bodies will either heat up or cool down toward 
the gas temperature of 150 K. The results are shown in 
the center panel of Figure 20, which illustrates Tg (top) 
and the normalized difference between numerical and an¬ 
alytical solutions (bottom) versus time. 

Finally, we discuss a test in which a particle is neither 
subject to frictional heating (as in the first test of this 
section) nor to radiative heating and cooling (as in the 
second test). We assume that a particle loses mass, due 
to ablation, releasing vaporization energy in the process. 
For simplicity, the variation of the particle mass is such 
that dRgjdt is constant, so that dMg/dt oc Thus, 
the energy budget reduces to 


Cs{Tg) 


d^ 

dt 



(B18) 


where, Lg , the specific energy of vaporization, is constant 
(see PPR88). The specific heat has the form Cg = cT^ 
{b ^ 1), which is an approximation to the specific heat of 
ice between 30 K and 273 K (Haynes 2011). The solution 
to Equation (B18) is 


Tg = 



Lg In 



i/(&+i) 


(B19) 


where Tgo is the initial particle temperature. For Tg > 
Tg, we have 


m) 


1 

4^3 
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In 


\Ts-tJ 


+ 2 arctan 



(B16) 


where T^q and Rgo are the initial temperature and radius 
of the particle. 

The right panel of Figure 20 shows the thermal evolu¬ 
tion of particles that ablate and lose energy. Again, it 
is assumed that the particle orbits the star with a pe¬ 
riod 27r/r^o- Tests are performed for various values T^o 
and i?so- The temperature is illustrated on top while the 
normalized difference between the computed and the an- 
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alytical temperature in Equation (B19) is shown on the 
bottom of the figure. 
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